# Improvements to the Search for Cosmic Dawn Using the Long Wavelength Array

Christopher DiLullo University of New Mexico
210 Yale Blvd. NE
Albuquerque, NM 87131, USA  Lab for Observational Cosmology, Code 665, NASA
Goddard Space Flight Center
8800 Greenbelt Rd.
Greenbelt, MD 20771, USA  Jayce Dowell University of New Mexico
210 Yale Blvd. NE
Albuquerque, NM 87131, USA  Gregory B. Taylor University of New Mexico
210 Yale Blvd. NE
Albuquerque, NM 87131, USA

###### Abstract

We present recent improvements to the search for the global Cosmic Dawn
signature using the Long Wavelength Array station located on the Sevilleta
National Wildlife Refuge in New Mexico, USA (LWA–SV). These improvements are
both in the methodology of the experiment and the hardware of the station. An
improved observing strategy along with more sophisticated temperature
calibration and foreground modelling schemes have led to improved residual RMS
limits. A large improvement over previous work using LWA–SV is the use of a
novel achromatic beamforming technique which has been developed for LWA–SV. We
present results from an observing campaign which contains 29 days of
observations between March $10^{\rm{th}}$, 2021 and April $10^{\rm{th}}$
2021\. The reported residual RMS limits are 6 times above the amplitude of the
potential signal reported by the Experiment to Detect the Global EoR Signature
(EDGES) collaboration.

Long Wavelength Array; Cosmic Dawn; First Light

##  1 Introduction

Detecting the formation of the first stars, known as Cosmic Dawn, and probing
the subsequent Epoch of Reionization (EoR) remain some of the biggest goals of
observational cosmology. The efforts to detect signals from these times in
cosmic history focus either on direct imaging of the first luminous sources in
the infrared, a major science goal of the upcoming James Webb Space Telescope
(Gardner et al., 2006), or detecting the redshifted 21 cm signal from neutral
hydrogen present during these periods (Furlanetto et al., 2006; Morales &
Wyithe, 2010; Pritchard & Loeb, 2012). The redshifted 21 cm signal offers the
only ground-based approach to studying Cosmic Dawn and the EoR since Earth’s
atmosphere is opaque in the relevant infrared band.

The 21 cm signal provides an especially strong probe which traces the
evolution of neutral hydrogen throughout the Epoch of Reionization. Measuring
the brightness of the signal across frequency would give insights into how
astrophysical conditions evolve as the first luminous sources begin to
dominate the Universe (Cohen et al., 2017). The 21 cm signal is generally
broken into two observable signals: a global signal which is detectable over
large angles on the sky and a spatially inhomogeneous signal which corresponds
to the growth of HII regions during the EoR. The global Cosmic Dawn signal is
expected to manifest as a small absorption feature in the average spectrum of
the sky below 100 MHz. Experiments to detect the 21 cm signal focus on
detecting either the global signal by measuring the average spectrum of the
sky and accurately modelling and removing foreground signals (Sokolowski et
al., 2015; Singh et al., 2018; Price et al., 2018; McKinley et al., 2020) or
the inhomogeneous signal by measuring the three dimensional k-space power
spectrum whose structure arises from the formation of ionized regions during
the Epoch of Reionization (Parsons et al., 2010; DeBoer et al., 2017; Mertens
et al., 2020; Yoshiura et al., 2021). Novel interferometric techniques have
also been developed for future use in detecting the global signal (McKinley et
al., 2012; Vedantham et al., 2015; McKinley et al., 2020).

The potential detection of a global absorption signal reported by the
Experiment to Detect the Global EoR Signature (EDGES) collaboration (Bowman et
al., 2018) has sparked much interest and debate due to its unexpected shape
and amplitude (Hills et al., 2018; Bradley et al., 2019; Singh & Subrahmanyan,
2019; Sims & Pober, 2019). If validated, the EDGES absorption signal would
imply that current cosmological models fail to accurately predict the
astrophysical processes which govern the 21 cm physics of neutral hydrogen
throughout Cosmic Dawn and the Epoch of Reionization. This hints to either new
physics in the form of interactions between dark matter and baryons (Barkana,
2018; Muñoz & Loeb, 2018; Berlin et al., 2018) or the presence of a previously
unaccounted for radio background (Feng & Holder, 2018; Mirocha & Furlanetto,
2018; Dowell & Taylor, 2018).

The implications of the unexpected nature of the absorption signal reported by
EDGES warrant independent validation using a different instrument. The Long
Wavelength Array (LWA; Taylor et al., 2012; Cranmer et al., 2017) located in
New Mexico, USA is a radio telescope whose frequency coverage allows it to
possibly validate the EDGES signal. Initial work using the LWA station located
on the Sevilleta National Wildlife Refuge, NM, USA (LWA–SV) has been
previously presented (DiLullo et al., 2020). In that work, a novel beamforming
approach was presented which allows for increased sensitivity compared to
single dipole experiments, improved foreground minimization, and in situ
astronomical calibration. It was concluded that the dependence of beam
shape111Throughout this work, the phrase “beam shape” is used to represent the
combination of both the forward gain and the angular structure of the main
lobe of the beam. with frequency is a major challenge and that a custom
beamforming framework would have to be developed for LWA–SV which would allow
for achromatic beamforming.

The work detailed here describes recent improvements to the methods described
in DiLullo et al. (2020) which have improved the residual RMS limits by an
average factor of 3. These include hardware improvements to LWA–SV, changes to
the observational strategy, improvements to the data analysis, and the
successful adoption of achromatic beamforming. The paper is structured as
follows: Section 2 details the recent changes to LWA–SV and to the
observational strategy; Section 3 details the achromatic beamforming
framework, discusses tests to compare the sensitivity of a standard LWA–SV
beam to that of an achromatic beam, and presents a software package which can
be used to simulate the beam pattern, both standard and custom, for a general
antenna array; Section 4 will present results from a recent observing campaign
which employs all the recent improvements; and Section 5 discusses these
results and looks towards future work.

##  2 Recent Improvements

The improvements to the methods described in DiLullo et al. (2020) fall into
two categories: methodological improvements to the observing strategy and data
analysis and improvements to the hardware and software of LWA–SV.

###  2.1 Observing Strategy and Data Analysis

The observing strategy laid out in DiLullo et al. (2020) focuses on the use of
simultaneous observations of a cold region on the sky (RA $9^{\rm h}\ 38^{\rm
m}\ 40.56^{\rm s}$, Dec $+30^{\degree}\ 49^{\prime}\ 1.4^{\prime\prime}$;
J2000), referred to as the Science Field, and Virgo A (RA $12^{\rm h}\ 30^{\rm
m}\ 49.42^{\rm s}$, Dec $+12^{\degree}\ 23^{\prime}\ 28^{\prime\prime}$;
J2000). Simultaneous observations of Virgo A allowed for astronomical
calibration which would help account for any time dependent variation in the
system. However, this did not account for the differences between the shapes
of the science and the calibration beams as they were pointed towards
different locations on the sky. In addition to the instantaneous differences
in beam shapes, the different declinations of the Science Field and Virgo A
meant that the science and calibration beams would trace different paths
across the sky. These two effects caused inaccuracies in the calibration of
the science beam.

A new Science Field center pointing (RA: $11^{\rm h}\ 00^{\rm m}\ 49.42^{\rm
s}$, Dec: $+12^{\degree}\ 23^{\prime}\ 28^{\prime\prime}$; J2000) has been
identified. The choice of coordinates still places the Science Field in a
large cold region on the sky, but guarantees that the science and calibration
beams trace out the same arc on the sky. The observations of the Science Field
and Virgo A are also no longer done simultaneously, but rather during times
when the targets align in terms of position on the sky. This ensures that the
derived calibration better accounts for the changing forward gain of the beam
throughout a single observation.

The data analysis methodology has been improved through the adoption of fully
time dependent temperature calibration and Bayesian foreground modelling. Time
dependent temperature calibration is achieved by modelling the beam pattern of
LWA–SV for every $5^{\textrm{th}}$ pointing throughout the run at every
$5^{\textrm{th}}$ observed frequency in the band and convolving the models
with the Global Sky Model (GSM; de Oliveira-Costa et al., 2008). This creates
a set of model spectra between the relevant Local Sidereal Time (LST) ranges
which can be linearly interpolated in both time and frequency at the same
temporal and spectral resolution of an observation to yield a model dynamic
spectrum of the Virgo A field which is used to calibrate the raw data. It
should be noted that the choice to use every $5^{\textrm{th}}$ pointing and
frequency is purely arbitrary and was chosen to reduce required compute time.
These models can be generated at higher temporal and spectral resolution
before interpolation, but this becomes quite computationally expensive and is
not expected to significantly change the results since the sky changes
smoothly over time and frequency and this smooth nature is accurately captured
through interpolation. An example of the model dynamic spectrum of the Virgo A
beam, interpolated at 40 ms time resolution and 19.1 kHz frequency resolution,
is shown in Figure 1.

Figure 1: Model dynamic spectrum of the Virgo A beam. These are generated by
simulating the beam pattern at every $5^{\textrm{th}}$ pointing for every
$5^{\textrm{th}}$ frequency and convolving these simulations with the Global
Sky Model at 7.5 s resolution to yield temperature spectra. These spectra are
then interpolated at the frequency and time resolution of an observation. The
median spectrum and drift are shown at the top and right of each plot,
respectively.

Bayesian foreground modelling allows for the posterior distributions of the
model parameters to be explored to search for correlations between foreground
model parameters and 21 cm model parameters. Correlations between the
foreground and 21 cm models would be non-physical and would imply issues with
our foreground models. We use the Python package emcee (Foreman-Mackey et al.,
2013) which uses the Affine Invariant Markov chain Monte Carlo (MCMC) Ensemble
Sampler (Goodman & Weare, 2010) to maximize the likelihood function, which
essentially measures the probability of a dataset given a certain set of model
parameters. Bayes’ theorem relates the posterior probability distribution,
$\mathcal{P}(\mathbf{\Theta}|\mathbf{D},\mathcal{H})$, of a set of parameters
$\mathbf{\Theta}$ given the data $\mathbf{D}$ and some model $\mathcal{H}$ to
the likelihood $\mathcal{L}(\mathbf{D|\mathbf{\Theta},\mathcal{H}})$ via

$\mathcal{P}(\mathbf{\Theta}|\mathbf{D},\mathcal{H})=\frac{\mathcal{L}(\mathbf{D|\mathbf{\Theta},\mathcal{H}})\Pi(\mathbf{\Theta}|\mathcal{H})}{\mathcal{Z}(\mathbf{D}|\mathcal{H})},$
(1)

where $\Pi(\mathbf{\Theta}|\mathcal{H})$ is the prior distribution that
encodes information about the model parameters and
$\mathcal{Z}(\mathbf{D}|\mathcal{H})$ is the evidence, which is the integral
of the likelihood function over the prior domain. We follow Harker et al.
(2012) and Bernardi et al. (2016) by assuming Gaussian measurement noise in a
single frequency channel which allows us to write the likelihood
$\mathcal{L}_{i}$ of observing the sky temperature $T_{sky}(\nu_{i})$ at
frequency $\nu_{i}$ as

$\mathcal{L}_{i}(T_{sky}(\nu_{i})|\mathbf{\Theta})=\frac{1}{\sqrt{2\pi\sigma^{2}(\nu_{i})}}e^{-\frac{(T_{sky}(\nu_{i})-T_{m}(\nu_{i},\mathbf{\Theta}))^{2}}{2\sigma^{2}(\nu_{i})}},$
(2)

where $T_{m}$ is the sky model and $\sigma$ is the standard deviation of the
instrumental noise is a given frequency channel given by

$\sigma(\nu_{i})=\frac{T_{sky}(\nu_{i})}{\sqrt{\Delta\nu\Delta t}},$ (3)

where $\Delta\nu$ and $\Delta t$ are the channel bandwidth and integration
time, respectively.

Given Equation 2, the log-likelihood of the entire temperature spectrum can be
represented as

$\rm{ln}\mathcal{L}(\mathbf{T}_{sky}|\mathbf{\Theta})=\sum_{i=1}^{N}\rm{ln}\mathcal{L}_{i}(T_{sky}(\nu_{i})|\mathbf{\Theta}),$
(4)

where we have switched to the log-likelihood due to computational efficiency.
The sky model is a superposition of the foreground contribution and the 21 cm
signal, written as

$T_{m}(\nu_{i})=T_{f}(\nu_{i})+T_{21cm}(\nu_{i}),$ (5)

but for the work presented in this paper we do not attempt to fit the 21 cm
signal as our residual root mean square (RMS) is still too large to
meaningfully fit a 21 cm signal (See Section 4). We use a slightly modified
foreground polynomial than that used in Bernardi et al. (2016); however, it is
still a $N^{\rm{th}}$ order log-polynomial given by

$T_{f}(\nu_{i})=\sum_{n=0}^{N}p_{n}\left[\textrm{log}_{10}\left(\frac{\nu_{i}}{\nu_{0}}\right)\right]^{n},$
(6)

where $\nu_{0}$ is a reference frequency that is set to the center of the
band. We choose a uniform prior for each model parameter, $p_{n}$, which
assigns equal probability to all real valued solutions.

###  2.2 Upgrades to LWA–SV

There have been a few upgrades to the hardware of LWA–SV which should help in
the search for Cosmic Dawn. The observational bandwidth of LWA–SV has been
improved to 20 MHz per tuning per beam. This new bandwidth allows for
uninterrupted frequency coverage in the range $50-85$ MHz with center tuning
choices of 60 and 75 MHz. Accounting for bandpass rolloff, we select the inner
80% of the band which yields a usable frequency range of $52-83$ MHz for the
output spectra. This frequency range should allow for the capture of the lower
edge of the EDGES signal through the center frequency of 78 MHz. This is a
large improvement over the spectra presented in DiLullo et al. (2020), which
had incomplete frequency coverage in the smaller range of $\approx 63-79$ MHz.

Along with the upgraded bandwidth capabilities of LWA–SV, a new weather
station was also deployed at the station which records the outside
temperature. We have previously seen that temperature variations inside the
electronics shelter can affect performance of the system, namely the analog
receiver (ARX) boards, and so we also expect outside temperature variations to
affect the performance of the front end electronics (FEE) on each dipole.
Outside temperature data from the new weather station can be combined with the
temperature data from the ARX boards to allow for a 2-dimensional fit which
better captures gain fluctuations than the previous simple 1-dimensional fit
with ARX temperature. However, we have found this effect to be so small in
magnitude that it is difficult to accurately remove this trend from the data
which is dominated by the changing sky.

The temperature change, both outside and in the electronics shelter, over a
single 1.5 hour observation is less than 1% compared to the change in sky
brightness due to the rising Galactic plane. The temperature dependence of the
FEE gain has been reported by Hicks et al. (2012) to be in the range 0.0042 –
0.0054 dB/°C between 20 – 100 MHz, which corresponds to a linear gain
difference of $\approx 0.001/\degree\rm{C}$. This small gain variation
combined with the small temperature changes over the course of a single
observation implies the FEE gain variation is negligible. The temperature
dependence of the ARX gain has not been directly measured, but we estimate it
to be of similar magnitude to the FEE temperature dependence.

Perhaps as we push our residuals below 1 K, this correction will become
necessary, but for now this approach does not seem to improve data quality.
Another improvement to note is the installation of second temperature sensor
for the HVAC controller within the electronics shelter. This has led to much
more stable temperatures within the shelter and hence more stable ARX
temperatures. We now see typical temperature variations on the order of
$\approx 2\degree$ F.

##  3 Achromatic Beamforming

A major challenge that all experiments attempting to detect Cosmic Dawn and
the EoR face is the frequency dependence of the antenna response. This
“chromaticity” of the receiving element produces spectral structure which can
obscure any potential cosmological signal. The array nature of LWA–SV helps as
beamforming allows for better rejection of foreground contributions by
focusing the main lobe in a region of the sky away from the Galactic plane;
however, this introduces a second chromaticity factor which can obscure the
signal. The beam pattern of an antenna array is not only a function of the
antenna gain pattern, which is chromatic, but also has intrinsic chromaticty.
This second chromaticity factor can be combated through a custom beamforming
framework which modifies the required complex beamforming coefficients in
order to tailor the array response to be more constant across the observed
frequency band. However, this still does not completely account for the
chromaticty of the gain pattern of the individual antennas within the array.

In DiLullo et al. (2020), a framework for custom beamforming at a single
frequency was presented. The beam pattern of LWA–SV can be shaped by modifying
the amplitude of the complex gains of each dipole in the array depending on
the frequency of interest and the desired pointing center on the sky. Model
custom beam patterns were presented and driftscan observations of Cygnus A (RA
$19^{\rm h}\ 59^{\rm m}\ 428.36^{\rm s}$, Dec $+40^{\degree}\ 44^{\prime}\
2.1^{\prime\prime}$; J2000) at transit confirmed that the beam full width at
half maximum (FWHM) could be predictably shaped using the framework. However,
the beam was only properly tuned for a single frequency. It was concluded that
in order to achieve near-achromatic beamforming, modifications would need to
be made to the beamformer pipelines of LWA–SV.

###  3.1 Implementation on LWA–SV

The beamformer pipelines are part of the digital backend of LWA–SV, known as
the Advanced Digital Processor (ADP; Cranmer et al., 2017; Price, 2017; Dowell
& Taylor, 2020). The ADP beamformer pipelines have been modified in order to
allow for pre-computed custom complex gains to be read in and overwrite the
complex gains which would normally be computed by the Monitor and Control
System at LWA–SV. Therefore, for an observation which consists of some number
of pointings which track a source on the sky, the necessary custom complex
gains can be computed for every observed frequency at each pointing. This
yields a set of custom complex gains which keep the beam pattern mostly
constant in both frequency and time as the beam tracks the source across the
sky. However, the beam pattern can still not be said to be completely constant
as the response pattern of the individual dipoles has both angular and
frequency dependence.

The major trade-off for the ability to shape the beam pattern of an antenna
array is a loss in sensitivity. The sensitivity of an array is given by its
System Equivalent Flux Density (SEFD)

$SEFD=\frac{k_{B}T_{sys}}{A_{e}}10^{26}\ \rm{Jy},$ (7)

where $k_{B}$ is the Boltsmann constant, $T_{sys}$ is the system temperature
in K, and $A_{e}$ is the effective area of the array in $\rm{m}^{2}$. The SEFD
is an estimate of the flux density of a source which would double the system
temperature, which means that a lower SEFD corresponds to a higher
sensitivity. $A_{e}$ is difficult to accurately measure for an antenna array,
but the simple estimate of

$A_{e}\approx N_{dip}\times A_{ant},$ (8)

where $N_{dip}$ is the number of antennas in the array and $A_{ant}$ is the
effective area of the antenna, is an idealized estimate. However, effects such
as mutual coupling of individual antennas can drastically change the effective
area of the array as a whole (Ellingson, 2007).

Compared to standard beamforming at LWA–SV where every antenna is equally
weighted, adjusting the amplitudes of the complex gains to shape the beam
effectively reduces the effective area of the entire array. Therefore, it is
expected for the SEFD of the array to increase as more dipoles are down
weighted to maintain beam shape with frequency. The SEFD of the array can be
estimated by observing a bright source, such as Cygnus A, and comparing the
observed power when the source is in and out of the beam. The SEFD can be
estimated with

$SEFD=\frac{S_{\nu}}{\frac{P_{on}}{P_{off}}-1},$ (9)

where $P_{on}$ is the observed power when the beam is pointed on source,
$P_{off}$ is the observed power when the beam is pointed off source, and
$S_{\nu}$ is the flux density of the source at the measured frequency, $\nu$.
This is typically done by pointing a beam at the transit position of a bright
source and letting the source drift into and out of the beam. Driftscans can
also be used to measure the FWHM of the beam main lobe. The major limitation
of driftscan observations is that only a single dimension of the beam can be
probed throughout the observation, so measuring the beam shape in both
dimensions on the sky requires either multiple beams or multiple days of
observation where the beam pointing center is below, on, and above the source
transit position so the source can probe different parts of the beam main
lobe.

A method to avoid the need for multiple days of observations is to
continuously repoint the beam in a basket weave pattern. This approach is
based on methods used by higher frequency arrays such as the Very Large Array.
The idea is to make two cuts along the source in different directions on the
sky. This allows for the 2-dimensional beam shape to be estimated at various
points in a single observation. However, at low frequencies ionospheric
scintillation on $\sim 10\ \rm{s}$ timescales can make things difficult since
the observed power at any position is a combination of the beam pattern and
the local ionospheric conditions. Thus, we modify the above approach to
interleave the off-source pointings with on-source pointings. These on-source
pointings serve as an ionospheric references which bracket each off-source
pointing and can be interpolated between to determine the amount of
ionospheric scintillation.

The modified basket weave approach is utilized on a weekly basis at LWA–SV in
order to monitor any changes in the station’s SEFD. These weekly observations
utilize standard beamforming on Cygnus A, so Cygnus A was a natural choice to
use when testing the sensitivity of the station with an achromatic beam. A
circular beam with FWHM of $4\degree$ was chosen for a basket weave
observation of Cygnus A; however, estimating the SEFD of LWA–SV turned out to
be challenging as the larger beam size meant that more power from the Galactic
plane was being picked up by the beam. This generally smoothed out the
contribution of Cygnus A to the total observed power and resulted in
unrealistic SEFD estimates. To avoid this, we used the beam-dipole observing
mode which LWA–SV offers to observe Cygnus A interferometrically using the
modified basket weave to determine the phase centers. Beam-dipole mode
correlates the beamformed output with the signal from a single antenna in the
array. We chose to correlate the beamformed output with the outrigger antenna
which is located $\sim 300\ \rm{m}$ to the West of the station center. This
two-element interferometer produces a fringe pattern on the sky which is a
function of the observed frequency, the baseline geometry, and the geometric
mean of the beam patterns of the station beam and the outrigger dipole. These
interferometric basket weave observations resolve out the diffuse structure
from the Galactic plane and help isolate the contribution from Cygnus A.
Results from observations using standard and achromatic beamforming can be
seen in Figure 2. The final SEFD estimates are averages between the two values
derived from the RA and Dec slices, respectively. The SEFD and FWHM estimates
for beam-dipole mode observations using standard and achromatic beamforming
are presented in Table 1. It is clear that the observations using achromatic
beams suffer a loss in sensitivity; however, the loss in sensitivity is less
than a factor of 2.

(a) Standard beamforming

(b) Achromatic beamforming

Figure 2: Interferometric basket weave results from observations of Cygnus A. The cuts along declination and right ascension are shown in each figure in the left and right columns, respectively. The lower and upper tunings are shown in each figure in the top and bottom rows, respectively. The achromatic beam shows less variation in FWHM across frequency in both directions on the sky. Table 1:  Beam-Dipole Mode Basket Weave Results Frequency | Standard Beamforming | Achromatic Beamforming
---|---|---
| SEFD: $4.0\ \rm{kJy}$ | SEFD: $4.8\ \rm{kJy}$
60 MHz |  |
| FWHM: $3^{\degree}22^{\prime}10.6^{\prime\prime}$ | FWHM: $5^{\degree}06^{\prime}09.3^{\prime\prime}$
| SEFD: $4.3\ \rm{kJy}$ | SEFD: $6.7\ \rm{kJy}$
75 MHz |  |
| FWHM: $2^{\degree}42^{\prime}16.1^{\prime\prime}$ | FWHM: $4^{\degree}44^{\prime}22.7^{\prime\prime}$

###  3.2 Beam Simulator: A Python Package to Simulate Array Beam Patterns

The achromatic beamforming work presented above might be of interest to a
broader community that is interested in either simply simulating the beam
pattern of a general antenna array or developing similar custom frameworks for
an antenna array. A Python package, aptly named Beam Simulator, has been
developed to help simulate the beam pattern of a general antenna array. It can
represent the elements of an antenna array down from a single cable up to the
entire array and is also compatible with output files from the Numerical
Electromagentics Code (NEC) which is used to numerically model the gain
pattern of antennas. This allows for the simulated beam pattern to encapsulate
the gain pattern of the individual array elements and not just account for the
geometry of the array. Beam Simulator also allows for custom beams to be
simulated using the same framework that has been developed for LWA–SV. There
are convenience functions built in to easily represent a LWA station, so this
should be of particular interest to any institutions that plan on hosting a
future LWA Swarm station (Taylor et al., 2019). However, it should also be of
interest to the more general community interested in representing the beam
pattern of any phased array. The package can be found and downloaded on GitHub
at https://github.com/cdilullo/beam_simulator.

##  4 Recent Observing Campaign and Results

An observing campaign was carried out between March $10^{\rm{th}},\ 2021$ and
April $10^{\rm{th}},\ 2021$ in order to build a sufficiently large data set. A
single observation consists of an achromatic beam with a main lobe FWHM of
$4{\degree}$ which tracks the Science Field (see Section 2.1) between the
local sidereal time (LST) range $10.28-11.78\ {\rm h}$ and a second achromatic
beam with a similar main lobe FWHM which tracks Virgo A between the LST range
$11.78-13.28\ {\rm h}$. We use the spectrometer observing mode offered by
LWA–SV which yields dynamic spectra with 40 ms time resolution and 1024
frequency channels each having 19.1 kHz of bandwidth. The individual
observations from each day are flagged for RFI in both time and frequency by
applying a median smoothing window to the median drift and spectrum,
respectively, to create smooth models of each and then flagging times or
frequencies which deviate from these models. This allows for the RFI flagging
to be controlled through 4 simple parameters: the smoothing window sizes and
the threshold values in both time and frequency. It was found that the
parameters which led to good RFI detection were a time smoothing window of
$30\ \rm{s}$, a frequency smoothing window of $250\ \rm{kHz}$, and time and
frequency thresholds of $3\textrm{-}\sigma$ and $5\textrm{-}\sigma$,
respectively, where $\sigma$ is the standard deviation of the data after the
smoothed model has been divided out. These parameters led to relatively
aggressive flagging for a single observation, $\approx 20\%$ of data for each
day, but this was intentional to try and capture low level RFI.

A common source of RFI in the relevant frequency band is ionospheric
reflections from terrestrial digital television (DTV) signals which manifest
as flares with a characteristic 6 MHz bandwidth. DTV channels $2\rm{-}6$ all
lie within our observed band and can be seen quite frequently. Due to the
short timescales of the bursts and the relatively small bandwidth that is
affected, these DTV flares are not captured in the above RFI flagging method
which is better suited at capturing long duration narrow-band RFI or short
duration wide-band RFI. To better capture these events, we search along the
pilot tone frequencies for each channel for times where the pilot tone is
$2\textrm{-}\sigma$ above the median value and flag the entire 6 MHz sub-band
corresponding to the DTV channel at these times.

After each day was flagged for RFI, we performed a moving median window
smoothing operation on the data in time in order to smooth out variations due
to ionospheric scintillation. Ionospheric scintillation causes random
variations in the brightness of our targets with time, which obscure the
temperature calibration. We chose a window size of 5 minutes in order to
properly smooth out scintillation. The individual flagged and smoothed
datasets were then temperature calibrated by multiplying a model temperature
dynamic spectrum of Virgo A (see Section 2.1) with the ratio between the raw
Science Field and Virgo A data. The calibrated Science Field data were then
integrated in time into 30 second bins before being combined into a single
dataset which contains the entire observing campaign. An extra RFI flagging
step was applied which flags any channel that has been flagged for more than
75% of the entire campaign. The fully flagged, calibrated, and integrated
dataset for the entire observing campaign is shown in Figure 3.

Figure 3: Calibrated datasets for the entire observing campaign at 30 second
integration time. The individual 1.5 hour observations can be seen as a
repeating pattern along the time axis. The full dataset consists of 29 days of
data taken between March $10^{\rm{th}},\ 2021$ and April $10^{\rm{th}},\
2021$.

We generated a measured spectrum by employing a bootstrapping algorithm to
generate samples from the full calibrated dataset. This bootstrapping
algorithm ensured that the output spectrum was the most “typical” and was not
influenced by random features, such as unflagged RFI, which were present in
only some of the integrated spectra. In order to generate the samples, we
first fit a smooth bandpass model to each integrated spectrum which was then
divided out to create a copy of the data which highlighted the unsmooth
features in each integrated spectrum. We then computed the standard deviation
across frequency for each of these residual spectra and computed the first and
third quartiles, $Q1$ and $Q3$, for the distribution of standard deviations.
We flagged outlier spectra whose standard deviation was above the standard
outlier threshold given by

$\sigma\geq 1.5\times(Q3-Q1).$ (10)

Out of a total of 5046 integrated spectra for the entire campaign, 405 outlier
spectra were identified for Tuning 1 XX, 274 for Tuning 1 YY, 163 for Tuning 2
XX, and 301 for Tuning 2 YY. The remaining non-outlier integrated spectra for
each tuning and polarization were then bootstrap sampled. The bootstrap
sampling randomly samples the set of non-outlier spectra and returns a set of
the same size. An average spectrum is then computed using the resampled set of
spectra. This was repeated 10,000 times in order to generate a set of average
spectra for each tuning and polarization. The average of these 10,000 spectra
was taken to be the “typical” spectrum. The two average spectra for a given
polarization were then combined to yield a single average spectrum which spans
the entire bandwidth of 52–83 MHz. A foreground model was then fit using the
Bayesian modelling framework described in Section 2.1. We chose to start 100
walkers around the zero vector and run a chain which was 25,000 steps long. It
was observed that the model struggled to properly fit the data across the
entire frequency band, with frequencies below $\approx 57.5$ MHz having the
poorest fit, so the lower edge of the band below 57.5 MHz was discarded before
the fitting. The average bootstrapped spectrum after trimming, along with the
fit foreground model and residuals, for each polarization is shown in Figure
4. Figure 5 shows a typical plot of the posterior distributions of the model
parameters for the XX polarization. A similar plot is generated for the YY
polarization.

Figure 4: Results of the MCMC fit. Figure 5: Posterior distributions for
foreground model parameters for Tuning 1 XX. The 1-dimensional posteriors show
the median value along with the $16^{\rm{th}}$ and $84^{\rm{th}}$ percentile
bounds. This plot was made with the Python package corner (Foreman-Mackey,
2016).

We also investigate the performance of maximally smooth functions (MSFs; Rao
et al., 2015, 2017) as foreground models. The Python package maxsmooth (Bevins
et al., 2021) is used to fit the same $5\textrm{th}$ order log-polynomial
foreground model (See Equation 6). We find that the MSF model performs well,
but there appears to be higher order spectral structure which is not captured.
This is not surprising as MSFs are designed to be maximally smooth so that
they preserve the non-smooth features of the data in residuals. The observed
spectrum along with the MSF model and residuals are shown in Figure 6. The
residual RMS across frequency for each polarization using both modelling
methods is summarized in Table 2.

Figure 6: Results of the maximally smooth function foreground model fit.  Table 2:  Residual RMS Across Frequency Foreground Modelling Method | XX | YY
---|---|---
MCMC | 2.47 K | 3.81 K
Maximally Smooth Function | 3.29 K | 5.26 K

##  5 Discussion and Conclusion

The residual RMS limits reported in Table 2 are encouraging as they are
beginning to approach the sub-Kelvin level that is required for a proper
validation of the reported EDGES signal. The desired goal is to attain a
residual RMS of $\approx 50\rm{mK}$ so that the EDGES signal could be
validated at the $10\textrm{-}\sigma$ level. However, pushing the residual RMS
to the sub-Kelvin level is not expected to be simple.

We believe that the largest systematic that is limiting our residual RMS is
the discrepancy between the beam model used in the temperature calibration and
the true beam pattern of LWA–SV on the sky. The model beam pattern predicts a
very smooth response which is unlikely to be realistic. We tested these
differences by taking an observation from a single day of the observing
campaign and fitting a smooth model to the data. Many spectral and temporal
structures remained after the smooth model was removed from the data which
implies that there are underlying low level non-smooth spectral structures
present in the data. These structures are most likely related to the sidelobe
pattern of the array which is difficult to measure and quantify.

There are also discrepancies between the model used for the gain pattern of a
single LWA dipole and the true gain pattern. The model is directly used in the
simulation of the beamformed pattern of the station as a whole and so is also
present in the temperature calibration. The model is based on Numerical
Electromagnetics Code (NEC) simulations of the LWA dipole, but directly
measuring it is extremely difficult at these low frequencies. There is an
ongoing collaboration with the External Calibrator for Hydrogen Observatories
(ECHO) team at Arizona State University whose goal is to measure the gain
pattern of the LWA dipole through use of a drone which transmits a tone of
known power at various altitudes and azimuths (Chang et al., 2015; Jacobs et
al., 2016).

Another concern is the presence of low level RFI which is difficult to detect
and flag via the median spectrum and drift method described in Section 4. Any
unflagged RFI will significantly increase the residual RMS, especially as we
approach the sub-Kelvin level. The bootstrapping algorithm described in
Section 4 should help reduce the contribution of random low level RFI which is
only present on some days of the observing campaign, but more systematic low
level RFI which is present throughout the entire observing campaign at the
same frequencies will still be problematic. As we decrease our residual RMS
limits, we will need to better understand the RFI environment of LWA–SV,
whether it is external or internal to the station itself.

Information about what physically might be limiting LWA–SV also can be gleaned
from the spectral structure of the residuals. There appears to be large scale
ripples in the residuals of size $\approx 20\rm{MHz}$ which could suggest some
physical signal reflection with a wavelength of 15 m, or even half this, is
present within the system. However, there is no characteristic length present
within the LWA–SV signal chain, so the origin of these ripples is still
unknown. Ripples such as these could possibly manifest from other effects such
as mutual coupling between antennas and cross coupling between signal paths on
the ARX boards, but these are much more difficult to isolate if they are
contributing to this structure.

It is also important to acknowledge that the above work does not calibrate out
additive noise from the electronics that is present in the data. This additive
noise will arise mainly from the amplifiers present in the FEE boards on each
antenna and the ARX boards which process the analog signal. Hicks et al.
(2012) has reported the noise temperature of the FEE boards to be in the range
255–273 K between 20–100 MHz, but the noise temperature of the ARX boards has
not been directly measured. The noise temperature of the ARX boards could be
measured by injecting a noise source of known temperature into the backend
electronics of LWA–SV. This was done in DiLullo et al. (2020), however the
purpose of that work was to verify that the system was stable over time and so
was not rigorously designed to measure the noise temperature of the ARX
boards. We leave a targeted study of the ARX noise temperature as a function
of frequency using signal injection for future work. The experience of other
experiments (Rogers & Bowman, 2012; Monsalve et al., 2017; Singh et al., 2018;
Nambissan et al., 2021) which have made similar measurements suggests this may
not be easy; however, our sky based calibration should capture the majority of
effects like this and so our noise measurements may not need to be so
painstaking.

The work presented above shows encouraging improvements to the previously
reported limits of LWA–SV. Future work will focus on better quantifying the
LWA dipole gain pattern and understanding discrepancies between models of the
station beam and the true beam pattern. This is crucial if the residual RMS is
to be reduced to the level where the 21 cm signal can be jointly fit with the
foreground model. A more robust RFI detection algorithm may have to be
developed in order to capture any low level RFI that is present throughout the
entire observing campaign.

Construction of the LWA has been supported by the Office of Naval Research
under Contract N00014-07-C-0147 and by the AFOSR. Support for operations and
continuing development of the LWA is provided by the Air Force Research
Laboratory and the National Science Foundation under grants AST-1835400 and
AGS-1708855. LWA

## References

  * Barkana (2018) Barkana, R. 2018, Nature, 555, 71
  * Berlin et al. (2018) Berlin, A., Hooper, D., Krnjaic, G., & McDermott, S. D. 2018, Physical review letters, 121, 011102
  * Bernardi et al. (2016) Bernardi, G., Zwart, J., Price, D., et al. 2016, Monthly Notices of the Royal Astronomical Society, 461, 2847
  * Bevins et al. (2021) Bevins, H., Handley, W., Fialkov, A., et al. 2021, Monthly Notices of the Royal Astronomical Society, 502, 4405
  * Bowman et al. (2018) Bowman, J. D., Rogers, A. E., Monsalve, R. A., Mozdzen, T. J., & Mahesh, N. 2018, Nature, 555, 67
  * Bradley et al. (2019) Bradley, R. F., Tauscher, K., Rapetti, D., & Burns, J. O. 2019, The Astrophysical Journal, 874, 153
  * Chang et al. (2015) Chang, C., Monstein, C., Refregier, A., et al. 2015, Publications of the Astronomical Society of the Pacific, 127, 1131
  * Cohen et al. (2017) Cohen, A., Fialkov, A., Barkana, R., & Lotem, M. 2017, Monthly Notices of the Royal Astronomical Society, 472, 1915
  * Cranmer et al. (2017) Cranmer, M. D., Barsdell, B. R., Price, D. C., et al. 2017, Journal of Astronomical Instrumentation, 6, 1750007
  * de Oliveira-Costa et al. (2008) de Oliveira-Costa, A., Tegmark, M., Gaensler, B., et al. 2008, Monthly Notices of the Royal Astronomical Society, 388, 247
  * DeBoer et al. (2017) DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, Publications of the Astronomical Society of the Pacific, 129, 045001
  * DiLullo et al. (2020) DiLullo, C., Taylor, G. B., & Dowell, J. 2020, Journal of Astronomical Instrumentation, 9, 2050008
  * Dowell & Taylor (2018) Dowell, J., & Taylor, G. B. 2018, The Astrophysical Journal Letters, 858, L9
  * Dowell & Taylor (2020) —. 2020, Long Wavelength Array Memo 214
  * Ellingson (2007) Ellingson, S. W. 2007, in 2007 IEEE Antennas and Propagation Society International Symposium, IEEE, 1501–1504
  * Feng & Holder (2018) Feng, C., & Holder, G. 2018, The Astrophysical Journal Letters, 858, L17
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306
  * Furlanetto et al. (2006) Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Physics reports, 433, 181
  * Gardner et al. (2006) Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Science Reviews, 123, 485
  * Goodman & Weare (2010) Goodman, J., & Weare, J. 2010, Communications in applied mathematics and computational science, 5, 65
  * Harker et al. (2012) Harker, G. J., Pritchard, J. R., Burns, J. O., & Bowman, J. D. 2012, Monthly Notices of the Royal Astronomical Society, 419, 1070
  * Hicks et al. (2012) Hicks, B. C., Paravastu-Dalal, N., Stewart, K. P., et al. 2012, Publications of the Astronomical Society of the Pacific, 124, 1090
  * Hills et al. (2018) Hills, R., Kulkarni, G., Meerburg, P. D., & Puchwein, E. 2018, Nature, 564, E32
  * Jacobs et al. (2016) Jacobs, D. C., Burba, J., Turner, L., & Stinnett, B. 2016, in 2016 IEEE Conference on Antenna Measurements & Applications (CAMA), IEEE, 1–3
  * McKinley et al. (2020) McKinley, B., Trott, C., Sokolowski, M., et al. 2020, Monthly Notices of the Royal Astronomical Society, 499, 52
  * McKinley et al. (2012) McKinley, B., Briggs, F., Kaplan, D. L., et al. 2012, The Astronomical Journal, 145, 23
  * Mertens et al. (2020) Mertens, F. G., Mevius, M., Koopmans, L. V., et al. 2020, Monthly Notices of the Royal Astronomical Society, 493, 1662
  * Mirocha & Furlanetto (2018) Mirocha, J., & Furlanetto, S. R. 2018, Monthly Notices of the Royal Astronomical Society, 483, 1980
  * Monsalve et al. (2017) Monsalve, R. A., Rogers, A. E., Bowman, J. D., & Mozdzen, T. J. 2017, The Astrophysical Journal, 835, 49
  * Morales & Wyithe (2010) Morales, M. F., & Wyithe, J. S. B. 2010, Annual review of astronomy and astrophysics, 48, 127
  * Muñoz & Loeb (2018) Muñoz, J. B., & Loeb, A. 2018, Nature, 557, 684
  * Nambissan et al. (2021) Nambissan, J., Subrahmanyan, R., Somashekar, R., et al. 2021, Experimental Astronomy, 1
  * Parsons et al. (2010) Parsons, A. R., Backer, D. C., Foster, G. S., et al. 2010, The Astronomical Journal, 139, 1468
  * Price et al. (2018) Price, D., Greenhill, L., Fialkov, A., et al. 2018, Monthly Notices of the Royal Astronomical Society, 478, 4193
  * Price (2017) Price, D. C. 2017, Long Wavelength Array Memo 213
  * Pritchard & Loeb (2012) Pritchard, J. R., & Loeb, A. 2012, Reports on Progress in Physics, 75, 086901
  * Rao et al. (2015) Rao, M. S., Subrahmanyan, R., Shankar, N. U., & Chluba, J. 2015, The Astrophysical Journal, 810, 3
  * Rao et al. (2017) —. 2017, The Astrophysical Journal, 840, 33
  * Rogers & Bowman (2012) Rogers, A. E., & Bowman, J. D. 2012, Radio Science, 47
  * Sims & Pober (2019) Sims, P. H., & Pober, J. C. 2019, arXiv preprint arXiv:1910.03165
  * Singh & Subrahmanyan (2019) Singh, S., & Subrahmanyan, R. 2019, arXiv preprint arXiv:1903.04540
  * Singh et al. (2018) Singh, S., Subrahmanyan, R., Shankar, N. U., et al. 2018, Experimental Astronomy, 45, 269
  * Sokolowski et al. (2015) Sokolowski, M., Tremblay, S. E., Wayth, R. B., et al. 2015, Publications of the Astronomical Society of Australia, 32
  * Taylor et al. (2012) Taylor, G., Ellingson, S., Kassim, N., et al. 2012, Journal of Astronomical Instrumentation, 1, 1250004
  * Taylor et al. (2019) Taylor, G., Dowell, J., Pihlström, Y., et al. 2019, Bulletin of the AAS, 51.  https://baas.aas.org/pub/2020n7i002
  * Vedantham et al. (2015) Vedantham, H., Koopmans, L., De Bruyn, A., et al. 2015, Monthly Notices of the Royal Astronomical Society, 450, 2291
  * Yoshiura et al. (2021) Yoshiura, S., Pindor, B., Line, J., et al. 2021, Monhtly Notices of the Royal Astronomical Society, 505, 4775


# Uncountably many cases of Filippov’s sewed focus

Paul Glendinning, S. John Hogan, Martin Homer,
Mike R. Jeffrey, Robert Szalai School of Mathematics, University of
Manchester, Manchester M13 9PL.Corresponding author:
<EMAIL_ADDRESS>of Engineering Mathematics, University of
Bristol, Bristol BS8 1TW.

###### Abstract

The sewed focus is one of the singularities of planar piecewise smooth
dynamical systems. Defined by Filippov in his book Differential Equations with
Discontinuous Righthand Sides (Kluwer, 1988), it consists of two invisible
tangencies either side of the switching manifold. In the case of analytic
focus-like behaviour, Filippov showed that the approach to the singularity is
in infinite time.

For the case of non-analytic focus-like behaviour, we show that the approach
to the singularity can be in finite time. For the non-analytic sewed centre-
focus, we show that there are uncountably many different topological types of
local dynamics, including cases with infinitely many stable periodic orbits,
and show how to create systems with periodic orbits intersecting any bounded
symmetric closed set.

##  1 Introduction

Piecewise smooth systems have a wide range of applications in mechanics,
engineering, control theory and biology [4, 10]. Many of the fundamental
results on these systems are collected in Filippov’s seminal book [7]. Some of
his results have been partly rediscovered (see [8] for a case study from
bifurcation theory). Other results are not well-known, or are themselves
incomplete.

A general planar piecewise smooth differential equation with switching
manifold $\Sigma:y=0$ is given by

$\begin{array}[]{lll}\dot{x}=P^{+}(x,y),&\dot{y}=Q^{+}(x,y),&\text{if
$y>0$},\\\ \dot{x}=P^{-}(x,y),&\dot{y}=Q^{-}(x,y),&\text{if
$y<0$},\end{array}$ (1)

where the functions $P^{\pm}$ and $Q^{\pm}$ can be smoothly extended to a
neighbourhood of $\Sigma$. We say that (1) is $C^{k}_{*}$ if $P^{\pm},Q^{\pm}$
are $C^{k}$, $k\in\mathbb{Z}^{+}\cup\\{\infty\\}$, and $C^{\omega}_{*}$ if
they are analytic.111Note that Filippov [7, p. 206] uses a more general
definition: systems are said to be $C^{k}_{*}$ if $P^{\pm}$ and $Q^{\pm}$ are
$C^{k}$ and if $\Sigma$ is $C^{k+1}$, fixed and the same for all systems in
the class.

Filippov [7, p. 220] describes a classification of singularities of (1) on
$\Sigma$. Assuming that the singularity is at the origin $\text{O}:=(0,0)$,
type 3 singularities occur if

$P^{\pm}(0,0)\neq 0,\quad Q^{+}(0,0)=Q^{-}(0,0)=0.$ (2)

One type 3 singularity is the sewed focus222 Also called a fused, merged or
stitched focus., which is made up of two invisible tangencies either side of
the switching manifold [7, p. 234]. Locally, trajectories in $y>0$ are
inverted parabolas and in $y<0$ they are parabolas, creating a flow that can
be focus-like (see Figure 1) or centre-like. The conditions for a sewed focus
are, up to a change in the direction of time,

$P^{+}(0,0)=b^{+}>0,\quad P^{-}(0,0)=b^{-}<0,$ (3)

and

$xQ^{\pm}(x,0)<0\quad{\rm if}~{}x\neq 0$ (4)

for $x$ in some neighbourhood $V$ of $x=0$. The simplest sufficient condition
to satisfy (4) locally is

$Q^{+}_{x}(0,0)<0,\quad Q^{-}_{x}(0,0)<0.$ (5)

Note that (4) and (5) each imply that $\dot{y}$ has the same sign at $(x,0)$
for all $x\neq 0$ on both sides of $y=0$. Hence (1) admits unique solutions
passing through $\Sigma$. There is no sliding333The vanishing of all
derivatives with respect to $x$ up to order $k-1$ and such that
${\partial^{k}Q^{\pm}}/{\partial x^{k}}<0$ also defines a sewed focus. and
hence no Filippov vector field.

\begin{overpic}[angle={90},height=202.3529pt]{Figure80.eps}
\put(106.0,44.8){\Large{$\Sigma:y=0$}} \put(37.0,44.7){\Huge{$\bullet$}}
\put(39.5,53.0){\Large{O}}
\put(33.7,35.3){{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\Huge{$\leftarrow$}}}
\put(43.0,58.8){{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\Huge{$\rightarrow$}}}
\put(5.0,52.0){\Large{$x_{1}$}} \put(11.2,45.7){\Large{$\bullet$}}
\put(86.0,52.0){\Large{$x_{2}$}} \put(85.0,45.7){\Large{$\bullet$}}
\put(65.0,38.0){\Large{$x_{0}$}} \put(66.0,45.5){\Large{$\bullet$}}
\end{overpic} Figure 1: Flow around a sewed focus at the origin O, formed from
two invisible tangencies, with switching manifold $\Sigma:y=0$, and
$x_{0,2}>0$, $x_{1}<0$. Reproduction of [7, Figure 80, p. 234].

If (1) is $C^{\omega}_{*}$, there are just two cases [7] up to time reversal:
the singularity is either i) a stable (unstable) focus, locally solutions tend
to $(0,0)$ as $t\to\infty$ ($t\to-\infty$), and solutions reach $(0,0)$ in
infinite time or ii) a centre and all solutions except the singularity are
periodic locally. Although these cases are given in [7], one aim of this paper
is to clarify the result and simplify the proof (see Theorem 1 below).

Another aim is to gain an understanding of the non-analytic cases, which were
not fully considered in [7]. In the case of focus-like behaviour we show that
_if the system is $C^{k}_{*}$, $k\in\mathbb{Z}^{+}\cup\\{\infty\\}$, then the
approach to the singularity can be in finite time_ (see Theorem 3 below).

We also consider more general sewed centre-focus444The centre-focus in smooth
systems is well-defined [12, p. 139]. behaviour, which was not considered at
all555“[T]he cases of [sewed] centre-focus are not considered here” [7, p.
222]. in [7]. In particular we show that there are uncountably many different
local structures666In fact, Filippov appears to have been aware of the
complexity possible in this case, stating that “…there are infinitely many
[sewed] [centre-foci]” [7, p. 250]. But he offers no proof of this statement.
The absence of any proof in the literature led us to results that follow. near
this singularity (see Theorem 6 below).

Our paper is organized as follows. In section 2 we review the theoretical
framework and results attributed by Filippov [7, p. 234] to Bautin and
Leontovich [3] and Skriabin [14], and give the description of the local
structure for piecewise analytic systems, $C^{\omega}_{*}$ (Theorem 1). This
theorem has been outlined before [7, §19.4]. It is included here because it
does not appear to be well known. Our proof clarifies and simplifies the one
given in [7, §19.4].

In section 3 we show that the approach to the singularity can be in finite
time for piecewise non-analytic systems, $C^{k}_{*}$,
$k\in\mathbb{Z}^{+}\cup\\{\infty\\}$. In section 4 we describe an example of a
sewed focus with infinitely many isolated stable periodic orbits (Theorem 4)
and show how this can be generalized to create systems with periodic orbits
intersecting any bounded symmetric closed set (Theorem 6).

##  2 The analytic case

If (1) is $C^{\omega}_{*}$, there are two types of local behaviour, up to time
reversal.

######  Theorem 1.

Suppose that the origin is a sewed focus of a $C^{\omega}_{*}$ system (1).
Then the origin is either i) a centre or ii) all solutions approach the
singularity in infinite time as either $t\to\infty$ (stable sewed focus) or
$t\to-\infty$ (unstable sewed focus).

###### Proof.

The approach of Bautin and Leontovich [3], as described by Filippov [7, p.
234], is to define functions $\sigma^{\pm}$ in a neighbourhood of the
singularity.

With the time convention as in (3), in $y\geqslant 0$, the function
$\sigma^{+}$ takes $(x_{1},0)$ on $\Sigma$, where $x_{1}<0$, to the next
intersection of the flow with $\Sigma$, given by $(x_{2},0)$, in forwards
time, where $x_{2}>0$. Hence in $x<0$, $\sigma^{+}(x_{1})=x_{2}$, locally; see
Figure 1. We can extend $\sigma^{+}$ to $x\geqslant 0$ by defining
$\sigma^{+}(0)=0$ and $\sigma^{+}(x_{2})=x_{1}$ with $x_{2}>0$ iff
$\sigma^{+}(x_{1})=x_{2}$.

The function $\sigma^{-}$ on $\Sigma$ following the flow in $y\leqslant 0$, is
defined in the same way: in $x<0$, $\sigma^{-}(x_{1})=x_{0}>0$ and in
$x\geqslant 0$, $\sigma^{-}(0)=0$ and $\sigma^{-}(x_{0})=x_{1}$.

Hence

$\sigma^{\pm}(\sigma^{\pm}(x))=x.$ (6)

We claim that $\sigma^{\pm}$ extend to analytic functions on a neighbourhood
of $x=0$ with

$\sigma^{\pm}(0)=0,\quad\sigma^{\pm\prime}(0)=-1.$ (7)

It is straightforward to show that $\sigma^{\pm\prime}(0)=-1$. Differentiate
both sides of (6) with respect to $x$, to get
$\sigma^{\pm\prime}(\sigma^{\pm}(x))\sigma^{\pm\prime}(x)=1$. By assumption
$\sigma^{\pm}(0)=0$. Hence $(\sigma^{\pm\prime}(0))^{2}=1$,
$\sigma^{\pm\prime}(0)=\pm 1$. Since $x\sigma^{\pm}(x)<0$, then
$\sigma^{\pm\prime}(0)=-1$.

We delay the justification of the claim to analyticity to the end of the
proof. Let

$\chi(x)\equiv\sigma^{+}(x)-\sigma^{-}(x)$ (8)

which is analytic on a neighbourhood of the origin. Using (7), we have

$\chi(0)=0,\quad\chi^{\prime}(0)=0.$ (9)

We begin by showing that if $\chi(x)>0$ ($<0$) $\forall x<0$ in a
neighbourhood of the singularity, then $\chi(x)>0$ ($<0$) $\forall x>0$
(locally) as well.

Clearly $\sigma^{\pm\prime}(x)<0$ on a neighbourhood of the origin so
$\sigma^{\pm}$ are decreasing functions777Also true since orbits of a vector
field on the plane can not intersect. of $x$ locally. Suppose that
$\chi(x)>0$, that is $\sigma^{-}(x)<\sigma^{+}(x)$, for all $x<0$ in a
neighbourhood of $x=0$. If we let $\sigma^{+}(x)=x_{+},\sigma^{-}(x)=x_{-}$,
then $x_{+}>x_{-}>0$. Since $\sigma^{+}$ is decreasing this implies that
$\sigma^{+}(x_{+})<\sigma^{+}(x_{-})$. Subtract $\sigma^{-}(x_{-})$ from both
sides and hence

$0=\sigma^{+}(x_{+})-\sigma^{-}(x_{-})<\sigma^{+}(x_{-})-\sigma^{-}(x_{-})$

where the equality holds because
$\sigma^{+}(x_{+})=\sigma^{-}(x_{-})=\sigma^{\pm}(\sigma^{\pm}(x))=x$ using
(6). Hence if $\chi(x)>0$ in $x<0$ locally then $\chi(x)>0$ in $x>0$ locally.
This follows from Figure 2, since the functions $z=\sigma^{\pm}(x)$ both
decrease with $x$ and are symmetric about the line $z=x$. If $\chi(x)<0$ then
the argument is completely analogous.

\begin{overpic}[height=177.06458pt]{Figure82.eps}
\put(95.0,25.0){$\sigma^{+}(x)$} \put(95.0,0.0){$\sigma^{-}(x)$}
\put(25.0,95.0){$\sigma^{+}(x)$} \put(-10.0,95.0){$\sigma^{-}(x)$}
\put(37.0,80.0){$z$} \put(85.0,37.0){$x$} \put(82.0,75.0){$z=x$} \end{overpic}
Figure 2: Demonstration that if $\chi(x)\equiv\sigma^{+}(x)-\sigma^{-}(x)>0$
$\forall x<0$ in a neighbourhood of the singularity, then $\chi(x)>0$ $\forall
x>0$ (locally) as well. Reproduction of [7, Figure 82, p. 235].

Since $\chi$, $\sigma^{+}$ and $\sigma^{-}$ are analytic then either i)
$\chi(x)\equiv 0$ on a neighbourhood of $x=0$, which implies that all local
solutions are periodic and so the origin is a centre, or ii) there exists
$r\geqslant 2$ such that $\chi^{(k)}(0)=0$, $k=1,\dots,r-1$ and
$\chi^{(r)}(0)\neq 0$. In this case the previous argument implies that $r$ is
even, $r=2p$ say, and so locally

$\chi(x)=\frac{\chi^{(2p)}(0)}{(2p)!}x^{2p}+o(x^{2p}).$ (10)

Hence in a neighbourhood of $x=0$ either $\chi(x)<0$ for all $x\neq 0$ (a
stable focus) or $\chi(x)>0$ for all $x\neq 0$ (an unstable focus).

The proof that the approach of the origin takes place in infinite time for the
sewed focus is given in [7, pp. 237–238]. Suppose that the first non-zero
derivative of $\chi$ at the origin is at order $2p$, as in (10). Let
successive intersections of a trajectory with the $x$-axis be denoted by
$x_{n}$. By definition (8)

$x_{2n+2}=x_{2n}+\chi(x_{2n+1}),$

where $x_{2n}>0,\,\forall\,n$. Then (7) implies that $x_{2n+1}/x_{2n}\to-1$ as
$n\to\infty$. Hence by (10)

$x_{2n+2}-x_{2n}=ax^{2p}_{2n+1}+o(x^{2p}_{2n+1})$

where $a={\chi^{(2p)}(0)}/{(2p)!}\neq 0$. Consider the stable888The unstable
sewed focus $a>0$ is similar after reversing the direction of time. sewed
focus ($a<0$). Then for sufficiently large $n$ there exists $A>0$ such that

$x_{2n+2}\geqslant x_{2n}-Ax_{2n}^{2p}\geqslant x_{2n}-Ax_{2n}^{2}.$

We will now show that $\sum_{1}^{\infty}x_{2n}$ diverges. Since the speed in
the $x$-direction is bounded away from zero in a neighbourhood of the origin,
then the time taken to approach the origin will also then diverge.

Let $u_{n+1}=u_{n}-Au_{n}^{2}$. If $u_{1}=x_{2}$ then
$\sum_{1}^{\infty}x_{2n}\geqslant\sum_{1}^{\infty}u_{n}$. Let
$v_{n}=\frac{x_{2}}{2n}$ for $n\geqslant 2$, chosen because
$\sum_{1}^{\infty}v_{n}$ diverges.

Now compare $u_{n}$ with $v_{n}$. Clearly $u_{1}>v_{1}$. Now suppose that
$u_{n}>v_{n}\,\forall\,n\in[2,m]$. Then

$u_{m+1}=u_{m}-Au_{m}^{2}>\frac{x_{2}}{2m}-A\frac{x_{2}^{2}}{4m^{2}}.$

So $u_{m+1}>v_{m+1}$ provided

$\frac{x_{2}}{2m}-A\frac{x_{2}^{2}}{4m^{2}}>\frac{x_{2}}{2(m+1)}$

or, on rearranging

$2m^{2}x_{2}>Ax_{2}^{2}m(m+1).$

Thus if $x_{2}$ is sufficiently small then $u_{m}>v_{m}$. Hence both
$\sum_{1}^{\infty}x_{2n}$ and $\sum_{1}^{\infty}u_{n}$ diverge, by the
comparison test.

It remains only to demonstrate that $\sigma^{\pm}$ are analytic on a
neighbourhood of $x=0$. Filippov [7] refers to classic analysis of Andronov et
al. [2], but here we sketch a rather simpler explanation due to Chillingworth
[5].

Consider $y\geqslant 0$ (the argument in $y<0$ is similar). We have, from (1),

$\dot{x}=P^{+}(x,y),\quad\dot{y}=Q^{+}(x,y)$

with, from (3) and (4),

$P^{+}(0,0)=b^{+}>0,\quad xQ^{+}(x,0)<0.$

Both $P^{+}$ and $Q^{+}$ are analytic in a neighbourhood of $(0,0)$. Then
write $P^{+}(x,y)=b^{+}+p^{+}(x,y)$ with $p^{+}$ analytic and $p^{+}(0,0)=0$.
Write $Q^{+}(x,y)=-x^{2k+1}r^{+}(x)+ys^{+}(x,y)$ for some $k\geqslant 0$, with
$r^{+}(0)=\frac{1}{(2k+1)!}\frac{\partial^{2k+1}Q^{+}}{\partial
x^{2k+1}}(0,0)\neq 0$ (an odd power is needed to ensure that $xQ^{+}(x,0)<0$
locally), with both $r^{+}$ and $s^{+}$ analytic near $(0,0)$.

Figure 3: Change of coordinates takes $y=0$ to a generalized parabola.

Then ${\mathrm{d}y}/{\mathrm{d}x}=F^{+}(x,y)$ where $F^{+}$ is locally
analytic and given by

$F^{+}(x,y)=\frac{Q^{+}(x,y)}{P^{+}(x,y)}=-x^{2k+1}R^{+}(x)+yS^{+}(x,y)$ (11)

with $R^{+}$ and $S^{+}$ analytic, close to $r^{+}$ and $s^{+}$ up to scaling.
Thus solutions must lie999Substitute into (11). on analytic curves of the form

$y=c-x^{2k+2}W^{+}(x)$

where $c$ is constant, $W^{+}$ is analytic, and $W^{+}(0)>0$.

Consider a flow box in a neighbourhood of the origin bounded by two
trajectories, one with $c>0$ and the other with $c<0$. By an analytic change
of coordinates $(x,y)\to(u,v)$ this can be transformed101010One example of a
simple transformation is $(x,y)\to(u,v)=(x,y+x^{2k+2}W^{+}(x))$. In that case,
$y=0$ would map to $v=W^{+}(u)u^{2k+2}$ and so $C(1+G(u))=W^{+}(u)$. to a flow
where solutions lie on horizontal lines as shown in Figure 3. Then the line
$y=0$ maps to a curve

$v=C(1+G(u))u^{2k+2}.$

An analytic near identity change of coordinates, leaving $u$ unchanged, of the
form $(u,v)\mapsto(u,w)$, brings this to

$w=u^{2k+2}.$

In these new coordinates the return map is simply $u\to-u$ which is analytic.
Hence $\sigma^{+}$ is analytic, concluding the proof.

∎

Provided the smallest non-zero derivative of $\chi$ exists, the proof also
works for $C^{k}$ functions.

######  Theorem 2.

If the origin is a sewed focus of a $C^{k}_{*}$ system, $k\geqslant 2$, and
for some $p$ such that $2p\leqslant k$ the function $\chi$ is $C^{2p}$ with

$\chi^{(s)}(0)=0,~{}s=1,2,\dots,2p-1,\quad\chi^{(2p)}(0)\neq 0,$

then locally trajectories approach the origin in infinite time, as
$t\to\infty$ if $\chi^{(2p)}(0)<0$ and as $t\to-\infty$ if $\chi^{(2p)}(0)>0$.

##  3 Finite time approach

The approach to the origin is asymptotic for $C^{k}_{*}$ systems satisfying
the conditions of Theorem 2. The aim of this section is to show that finite
time approach is also possible for these non-analytic systems.

######  Theorem 3.

For all $k\in\mathbb{Z}^{+}\cup\\{\infty\\}$, there exists a $C^{k}_{*}$
system of the sewed focus and a neighbourhood $U$ of the origin such that if
$(x_{0},y_{0})\in U$ then the trajectory through $(x_{0},y_{0})$ tends to
$(0,0)$ in finite time.

###### Proof.

First we assume $k\in\mathbb{Z}^{+}$. Let $H(x)$ denote the Heaviside
function,

$H(x)=\begin{cases}1&{\rm if}~{}x\geqslant 0\\\ 0&{\rm if}~{}x<0.\end{cases}$
(12)

Consider system (1) with

$\displaystyle P^{+}(x,y)$ $\displaystyle=1$ $\displaystyle Q^{+}(x,y)$
$\displaystyle=-4x^{3}H(x)-8x^{7}H(-x),$ (13) $\displaystyle P^{-}(x,y)$
$\displaystyle=-1,$ $\displaystyle Q^{-}(x,y)$
$\displaystyle=-4x^{3}H(-x)-8x^{7}H(x).$

The symmetry and the values of the constants (other than their signs) are not
important, and have been chosen for algebraic convenience. This example is
$C^{2}_{*}$ and has a stable sewed focus at the origin. There is no sliding.

Given an initial point $(-x_{0},0)$ with $x_{0}>0$ sufficiently small, the
solution will strike $\Sigma$ at points $(x_{1},0),(-x_{2},0),\dots$ with
$x_{i}>0$. The total time $T$ taken is given by

$T=x_{0}+2\sum_{r=1}^{\infty}x_{r},$

since from (13) the speed in the $x$-direction is equal to 1.

Thus if $T<\infty$ the singularity is contained in the approaching trajectory,
and is reached in finite time. The trajectories of systems in the same
topological class are mapped to each other. So the inclusion of the origin O
in an approaching trajectory implies the existence of a new topological class,
not possible in the analytic case.

If $y>0$ then integral curves are given by

$y=c^{-}-x^{8},\ \ x<0,\quad y=c^{+}-x^{4},\ \ x>0$ (14)

where the constants $c^{\pm}$ are determined by initial conditions. The
integral curves in $y<0$ are determined by symmetry.

Thus a solution starting at $(-x_{0},0)$, $x_{0}>0$, lies on the integral
curve

$y=x_{0}^{8}-x^{8}$

and strikes the $y$-axis after time $x_{0}$ at $y=x_{0}^{8}$. It now continues
on an integral curve of the form $y=c^{+}-x^{4}$. Since $y=x_{0}^{8}$ when
$x=0$, $c^{+}=x_{0}^{8}$. This integral curve in turn strikes $\Sigma$ at
$(x_{1},0)$, $x_{1}>0$, where $x_{0}^{8}-x_{1}^{4}=0$, i.e. $x_{1}=x_{0}^{2}$,
after a further $x_{1}$ units of time.

Now in $y<0$, which by symmetry is effectively equivalent, the solution
through $(x_{1},0)$ strikes $\Sigma$ at $(-x_{2},0)$ where
$x_{2}=x_{1}^{2}=x_{0}^{4}$ after time $x_{1}+x_{2}$. By induction, the
infinite sequence of intersections is given by the points
$((-1)^{(r+1)}x_{r},0)$ where

$x_{r}=x_{0}^{2^{r}}$

and hence

$T=x_{0}+2\sum_{r=1}^{\infty}x_{0}^{2^{r}}$ (15)

The sum (15) is certainly less than twice the sum of all powers of $x_{0}$,
which converges if $x_{0}<1$ as it is a geometric progression, so this sum
also converges to a finite value.

Replacing $4x^{3}$ and $8x^{7}$ in (13) by $2kx^{2k-1}$ and $4kx^{4k-1}$,
respectively, with $k\geqslant 3$, the relation between subsequent
intersections is unchanged and so this example can be made $C^{r}_{*}$ for any
$r<\infty$.

We now give an example111111We are extremely grateful to an anonymous reviewer
for this example. of a $C^{\infty}_{*}$ system with finite time approach to
the origin. Consider system (1) with

$\displaystyle P^{+}(x,y)$ $\displaystyle=1,$ $\displaystyle Q^{+}(x,y)$
$\displaystyle=\begin{cases}-\frac{2}{x^{3}}e^{-1/x^{2}},&x<0,\\\ 0,&x=0,\\\
-\frac{1}{x^{2}}e^{-1/x},&x<0\end{cases}$ (16) $\displaystyle P^{-}(x,y)$
$\displaystyle=-1,$ $\displaystyle Q^{-}(x,y)$
$\displaystyle=\begin{cases}\frac{1}{x^{2}}e^{-1/x},&x<0\\\ 0,&x=0,\\\
-\frac{2}{x^{3}}e^{-1/x^{2}},&x>0.\end{cases}$

Then $x_{r}=x^{2}_{r-1}$ for all $r\geqslant 1$ and $T$, also given by
equation (15), is finite.

∎

Hence there is a new topological class with finite time approach to the
singularity for piecewise non-analytic systems, $C^{k}_{*}$ with
$k\in\mathbb{Z}^{+}\cup\\{\infty\\}$.

##  4 The general sewed focus

In this section, we will prove that there are three possibilities for the
local dynamics near a sewed focus, independent of the time taken to approach
the singularity121212The issue of finite versus infinite time approach in case
(a) is not included here as it was considered in section 3.. Also we show how
to construct a non-analytic piecewise smooth system with infinitely many
isolated stable periodic orbits.

######  Theorem 4.

Consider an isolated sewed focus singularity of a $C^{k}_{*}$ system (1),
where either $1\leqslant k\leqslant\infty$ or $k=\omega$. Let $\chi$ be
defined by (8). Then the singularity is either

  1. (a)

a stable (or an unstable) sewed focus (with $\chi(x)\neq 0$ for all $x\neq 0$
locally) or

  2. (b)

a sewed centre (with $\chi(x)\equiv 0$ locally) or

  3. (c)

a sewed centre-focus (where $\exists~{}x\neq 0$ arbitrarily close to $x=0$
such that $\chi(x)=0$).

If $k=\omega$ then only cases (a) or (b) can occur, otherwise all three cases
can occur.

###### Proof.

For $C^{\omega}_{*}$ systems, by the results of section 2, either i)
$\chi(x)\equiv 0$ (a sewed centre) or ii) $\exists~{}k\geqslant 1$ such that
$\chi^{(2k)}(0)\neq 0$ and all lower derivatives are zero. In the latter case,
either $\chi(x)>0$ or $\chi(x)<0$ for all $x\in B(0,\delta)\backslash\\{0\\}$
for $\delta>0$ small enough, giving either an unstable or a stable sewed focus
respectively. Nothing else is possible.

For $C^{k}_{*}$ systems, $1\leqslant k\leqslant\infty$ we consider vector
fields in $y<0$ of the form

$P^{-}(x,y)=-1,\quad Q^{-}(x,y)=-2x.$ (17)

In $y<0$ solutions lie on parabolas of the form $y=-c^{2}+x^{2}$. In $y>0$
define

$P^{+}(x,y)=1,\quad
Q^{+}(x,y)=-2x+2kx^{2k-1}f(\textstyle{\frac{1}{x}})-x^{2k-2}f^{\prime}(\textstyle{\frac{1}{x}})$
(18)

with $f(x)$ odd. By design, $Q^{+}$ is $C^{k}_{*}$ for sufficiently
differentiable $f$. So (1), with (17) and (18) is $C^{k}_{*}$ and satisfies
(2).

The flow in $y>0$ has been chosen so that solutions through $(-x_{0},0)$ lie
on curves

$y(x)=-x^{2}+x^{2k}f(\textstyle{\frac{1}{x}})+x_{0}^{2}+x_{0}^{2k}f(\textstyle{\frac{1}{x_{0}}})$
(19)

where we use the fact that $f$ is odd. Hence

$y(x_{0})=2x_{0}^{2k}f(\textstyle{\frac{1}{x_{0}}}).$

The zeros of this function correspond precisely to the closed orbits of the
vector field, provided that $y(x)>0$ for all $x\in(-x_{0},x_{0})$. Thus if we
choose

$f(x)=\sin 2\pi x$

and $k\geqslant 2$, then $\forall\,n\in\mathbb{R}$, $n\geqslant 4$, we have

$x_{0}=\frac{1}{n},$

which lies on a (different) periodic solution for each $n$; those with $n$
even are asymptotically stable and those with $n$ odd are unstable.

Choosing different forms of $f$ gives different zero sets. Hence case $(c)$
comes in infinitely many different types, depending on whether the components
of sets with $\chi(x)=0$ are intervals or points, and how these sets are
connected. ∎

We can now make some generalizations. We begin by recalling the following
result.

######  Lemma 5.

Let $E$ be any symmetric compact subset of $\mathbb{R}$ not containing the
origin (i.e. $E=-E$ and $0\notin E$), and
$k\in\mathbb{Z}^{+}\cup\\{\infty\\}$. Then there exists an odd $C^{k+1}$
function $f_{E}$ which, if $k<\infty$, is not a $C^{k+2}$ function such that
$f_{E}(x)=0$ if and only if $x\in E\cup\\{0\\}$.

###### Proof.

Given a symmetric compact set $E$ and $0\notin E$, let $E_{0}=E\cup\\{0\\}$
and $d(x,E_{0})=\min_{y\in E_{0}}\lvert x-y\rvert$. Define

$g_{E}(x)=-\textrm{sign}(x)[d(x,E_{0})]^{k+2}.$ (20)

Clearly $g_{E}$ is a continuous, odd function. Moreover, $g_{E}(x)=0$ iff
$x\in E_{0}$. Since $E_{0}$ is closed its complement is open and hence is a
countable union of disjoint open intervals. The function $g_{E}$ is $C^{k+1}$
except at the mid-points of each of these intervals. But $g_{E}$ can be
modified in a neighbourhood of each of these mid-points so that the
modification $f_{E}$ is $C^{k+1}$ there and still odd. This proves the first
part of the result if $k<\infty$.

If $k=\infty$ it is still the case that the complement of $E_{0}$ is open and
is thus a union of disjoint open intervals. This union is countable (as each
open interval contains a rational) and so we can label the intervals
$J_{i}=(a_{i}-\delta_{i},a_{i}+\delta_{i})$, so that the centre points are
$a_{i}$ and the length of $J_{i}$ is $2\delta_{i}$, $i=1,2,\dots$. Now let
$\psi(x)$ denote the standard $C^{\infty}$ bump function on the unit interval:

$\psi(x)=\begin{cases}\exp\left(-\frac{1}{1-x^{2}}\right)&\text{if $\lvert
x\rvert<1$}\\\ 0&\text{if $\lvert x\rvert\geqslant 1$},\end{cases}$ (21)

and define

$f_{E}(x)=\sum_{i=1}^{\infty}-\textrm{sign}(a_{i})\delta_{i}\psi\left(\frac{x-a_{i}}{\delta_{i}}\right).$
(22)

Then $f_{E}$ is $C^{\infty}$ and zero if and only if $x\in E_{0}$ ($0\in
E_{0}$ by definition so $a_{i}\neq 0$). Finally, $f_{E}$ is an odd function,
since $E_{0}$ is symmetric.

∎

The separation of $x=0$ from $E$ in the statement of Lemma 5 is because the
interpretation of the points in $E$ for the sewed focus (they are on periodic
orbits) is different from the sewed focus itself which has $x=0$.

Hence we have the rather remarkable theorem.

######  Theorem 6.

Let $E$ be any symmetric compact subset of $\mathbb{R}$ not containing the
origin. Then there exists a $C^{k}_{*}$ Filippov system (1), $1\leqslant
k\leqslant\infty$, with an isolated type 3 singularity, for which $(x_{0},0)$
lies on a periodic orbit if and only if $x_{0}\in E$.

_Proof:_ Choose $1\leqslant k\leqslant\infty$ and let $f_{E}(x)$ be the odd
$C^{k+1}$ or $C^{\infty}$ function of Lemma 5. Define the $C^{k}_{*}$ Filippov
system using (17) in $y<0$ and

$P^{+}(x,y)=1,\quad Q^{+}(x,y)=-2x+x^{2}f_{E}^{\prime}(x)+2xf_{E}(x)$ (23)

in $y>0$, which satisfies (2) since $f_{E}^{\prime}$ is continuous.
Furthermore, as $f_{E}^{\prime\prime}$ is continuous and $f_{E}(0)=0$,
$Q^{\pm}_{x}(x,y)<0$ in a neighbourhood of the origin. Hence in this
neighbourhood $Q^{\pm}(x,y)>0$ for $x<0$ and $Q^{\pm}(x,y)<0$ for $x>0$, which
guarantees that there is no sliding in the neighbourhood of the origin. Thus
solutions in $y>0$ through $(-x_{0},0)$ lie on sets

$y(x)=-x^{2}+x^{2}f_{E}(x)+x_{0}^{2}+x_{0}^{2}f_{E}(x_{0})$ (24)

by integrating ${\mathrm{d}y}/{\mathrm{d}x}={Q^{+}}/{P^{+}}$, using (23).

The flow in $y<0$ implies that an orbit through $(x_{0},0)$, $x_{0}>0$, next
intersects the switching manifold at $(-x_{0},0)$ and hence a periodic orbit
through $(\pm x_{0},0)$ exists iff the flow in $y>0$ takes $(-x_{0},0)$ to
$(x_{0},0)$, i.e. that $\sigma^{+}(-x_{0})=x_{0}$.

First suppose that $x_{0}>0$ and $(x_{0},0)$ lies on a periodic orbit. Then
following the flow in $y<0$, $(-x_{0},0)$ lies on the same periodic orbit, and
hence $y(\pm x_{0})=0$. Then (24) gives

$y(x_{0})=2x_{0}^{2}f_{E}(x_{0})=0,$

and so $x_{0}\in E$.

To prove sufficiency, suppose that $x_{0}>0$ and $x_{0}\in E$. Since $E$ is
symmetric $-x_{0}\in E$ too, and so $f_{E}(\pm x_{0})=0$. Then, recalling that
$f_{E}$ is odd, (24) gives

$y(x_{0})=2x_{0}^{2}f_{E}(x_{0})=0,$

so an orbit through $(-x_{0},0)$ crosses the switching manifold at
$(x_{0},0)$. Suppose this is not the first crossing, i.e. that there exists at
least one root $y(x^{\star})=0$ with $x^{\star}\in(-x_{0},x_{0})$. For all
$x\in E_{0}$, $\lvert x\rvert<x_{0}$,

$y(x)=-x^{2}+x_{0}^{2}+x^{2}f_{E}(x)=-x^{2}+x_{0}^{2}>0,$ (25)

so $x^{\star}$ must be in the complement of $E_{0}$, in one of the intervals
$J_{i}$, see Lemma 5. Equation (25) also implies that $y(x)>0$ at the
endpoints of all these intervals; since $y(x)$ is continuous, and has no
tangencies (as $Q^{\pm}(x,y)\neq 0$ for $x\neq 0$), there must be an even
number of roots of $y(x)$ in the interval $J_{i}$ containing $y^{\star}$. But
this contradicts the fact that $Q^{\pm}(x,y)$ do not change sign on any
interval $J_{i}$, as noted above. So there is no root
$x^{\star}\in(-x_{0},x_{0})$, and the trajectory in $y>0$ through $(-x_{0},0)$
next crosses the switching manifold at $(x_{0},0)$, and hence is a periodic
orbit for any $x_{0}\in E$.

Finally, let $x_{0}>0$ with $x_{0}\notin E$. Again, by symmetry of $E$,
$-x_{0}\notin E$ and this time (24) gives

$y(x_{0})=2x_{0}^{2}f_{E}(x_{0})\neq 0,$

so the trajectory through $(-x_{0},0)$ cannot be a periodic orbit.

∎

In Theorem 4, we have shown that there are uncountably many topological
classes of the sewed centre-focus, a type of non-analytic isolated sewed focus
singularity of a $C^{k}_{*}$ system (1), where $1\leqslant k\leqslant\infty$.
This result was generalized in Theorem 6.

##  5 Conclusion

In this paper we have analyzed the sewed focus, a singularity of planar
piecewise smooth systems. Although some analysis of this case is presented by
Filippov [7], we believe that many of our results are new. In the case of non-
analytic focus-like behaviour we show that the approach to the singularity can
be in finite time, forwards or backwards as appropriate. More generally, for
the non-analytic sewed centre-focus, we show that there are uncountably many
different topological types of local dynamics, including cases with infinitely
many stable periodic orbits.

The proof of our results uses a simple construction. By considering particular
cases in greater detail it is then possible to define $C^{k}_{*}$ systems,
$k\in\mathbb{Z}^{+}\cup\\{\infty\\}$, which have _infinitely many stable
periodic orbits_. This is another case of ‘infinitely many sinks’ arising in
piecewise smooth systems, although the mechanism is totally different from the
examples for the border collision normal form [13]. See also recent work [9],
which shows that the number of limit cycles for regularized piecewise
polynomial systems is unbounded.

Hence Hilbert’s 16th Problem, concerning the number of stable periodic orbits
of planar differential equations, is not relevant for the $C^{k}_{*}$ sewed
focus. This problem has attracted considerable recent attention in the case of
piecewise linear systems with a linear switching manifold [6, 11]. In
contrast, our paper answers the question for more general piecewise smooth
systems.

The sewed focus is one of the earliest examples of type 3 singularities to be
considered in the applied literature, with examples from the 1930s being
included in Andronov et al. [1, p. 511]. These examples are $C^{\omega}_{*}$,
so our more interesting results described here do not apply. But the effects
described may arise in examples with more complicated switches which change
higher derivatives of the flow in either $y>0$ or $y<0$ or both, as shown in
section 3.

_Acknowledgements:_ We are grateful to David Chillingworth (University of
Southampton) for his insights on the analyticity of return maps, and to James
Montaldi and Jonathan Fraser (University of Manchester) for helpful
discussions, especially concerning Lemma 5. We are extremely grateful to an
anonymous reviewer who supplied us with example (16), to complete the proof of
Theorem 3. SJH would like to thank the Hungarian Academy of Sciences for the
award of a Distinguished Guest Scientist Fellowship, during which part of this
paper was written.

## References

  * [1] A.A. Andronov, A.A. Vitt and S.E. Khaikin (1966) _Theory of Oscillators_ , Pergamon Press, Oxford.
  * [2] A.A. Andronov, E.A. Leontovich, I.I. Gordon and A.G. Maier (1971) _The Theory of Bifurcations of Dynamic Systems in the Plane_ , Israel Program for Scientific Translations, Jerusalem.
  * [3] N.N. Bautin and E.A. Leontovich (1976) _Methods and Techniques of Qualitative Study of Dynamical Systems in the Plane_ , Nauka, Moscow (in Russian).
  * [4] M. di Bernardo, C. Budd, A.R. Champneys and P. Kowalczyk (2008) _Piecewise-smooth Dynamical Systems Theory and Applications_ , Applied Mathematical Sciences, Vol. 163, Springer, London.
  * [5] D. Chillingworth (2016) private communication.
  * [6] R.D. Euzébio and J. Llibre (2015) On the number of limit cycles in discontinuous piecewise linear differential systems with two pieces separated by a straight line, _J. Math. Anal. Appl._ , 424, 475–486.
  * [7] A. F. Filippov (1988) _Differential Equations with Discontinuous Righthand Sides_ , Kluwer Academic Publishers, Dordrecht.
  * [8] S.J. Hogan, M.E. Homer, M.R. Jeffrey and R. Szalai (2016) Piecewise smooth dynamical systems: the case of the missing boundary equilibrium bifurcations, _J. Nonlinear Sci._ , 26, 1161–1173.
  * [9] R. Huzak and K.U. Kristiansen (2022) The number of limit cycles for regularized piecewise polynomial systems is unbounded, _J. Diff. Eqns_ , 342, 34–62.
  * [10] M.R. Jeffrey (2018) _Hidden Dynamics. The Mathematics of Switches, Decisions and Other Discontinuous Behaviour_ , Springer Nature, Switzerland AG.
  * [11] J. Llibre and X. Zhang (2018) Limit cycles for discontinuous planar piecewise linear differential systems separated by one straight line and having a center, _J. Math. Anal. Appl._ , 467, 537–549.
  * [12] L. Perko (2001) _Differential Equations and Dynamical Systems. 3rd Ed_ , Texts in Applied Mathematics, Vol. 7, Springer, New York.
  * [13] D.J.W. Simpson (2014) Sequences of periodic solutions and infinitely many coexisting attractors in the border-collision normal form, _Int. J. Bif. & Chaos_ 24:1430018.
  * [14] B.N. Skriabin (1978) The formation of a limit cycle at a “merged focus”, _PMM J. Appl. Math. Mech._ , 42, 1042–1045.