# CSS1603+19: a low-mass polar near the cataclysmic variable period minimum

Yiqi Liu,1 Hsiang-Chih Hwang,2 Nadia L. Zakamska1,2 and John R. Thorstensen 3
1Department of Physics & Astronomy, Johns Hopkins University, 3400 N Charles
St, Baltimore, MD 21218, USA
2School of Natural Sciences, Institute for Advanced Study, Princeton, 1
Einstein Drive, NJ 08540, USA
3Department of Physics and Astronomy, Dartmouth College, Hanover, NH 03755,
USA Email<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

CSS1603+19 is a cataclysmic variable (CV) with an orbital period of 81.96 min,
near the minimal period of cataclysmic variables. It is unusual in having a
strong mid-infrared excess inconsistent with thermal emission from a brown
dwarf companion. Here we present time-resolved multi-wavelength observations
of this system. WISE photometry indicates that the mid-infrared excess
displays a one-magnitude eclipsing-like variability during the orbit. We
obtained near-infrared and optical spectroscopy using Gemini, MDM and APO
telescopes. Near-infrared spectra show possible cyclotron features indicating
that the white dwarf has a magnetic field of about 5MG. Optical and near-
infrared spectra display double-peaked emission lines, with both components
showing strong radial velocity variations during the orbital period and with
the broad component leading the narrow component stably by about 0.2 of the
orbital phase. We construct a physical model informed by existing observations
of the system and determine that one component likely originates from the
accretion column onto the magnetized white dwarf in synchronous rotation with
the orbital motion and the other from the Roche overflow point. This allows us
to constrain the masses of the binary components to be $M_{1}>0.24M_{\odot}$
for the white dwarf accretor and $M_{2}=0.0644\pm 0.0074M_{\odot}$ for the
donor. We classify the system as an AM Herculis star, or a polar. It has
likely completed its stint on the period gap, but has not yet gone through the
period bounce.

###### keywords:

binaries: close – binaries: spectroscopic – novae, cataclysmic variables –
stars: low mass – white dwarfs

††pagerange: CSS1603+19: a low-mass polar near the cataclysmic variable period
minimum–References

##  1 Introduction

Cataclysmic variables (CVs) are compact binaries consisting of an accreting
white dwarf and a Roche-lobe filling companion star, or ‘donor’ (Warner,
1995). The formation of CVs is driven by angular momentum dissipation of a
post-common-envelop binary, where the orbital period of the system gradually
decreases from an initial orbital period P${}_{\rm orb}\gtrsim 10$ hours to an
orbital period minimum of P${}_{\rm min}\simeq 80$ min (Patterson, 1984;
Ritter & Kolb, 2003; Knigge et al., 2011; Kalomeni et al., 2016).

The complex distribution of observed orbital periods of CVs is accounted for
by the disrupted magnetic braking scenario (Rappaport et al., 1983). At
$P_{\rm orb}\gtrsim 3$ hours, the angular momentum loss for CVs is dominated
by magnetic braking (Patterson, 1984). In this scenario, the low-mass donor’s
magnetic field drives ionized stellar wind, which co-rotates with the star’s
field out to the Alfvénic limit, exerting a resistant torque on the donor star
and dissipating the angular momentum. In the meanwhile, in CVs the donor star
is tidally locked and is co-rotating with the system’s orbit. Thus, the
angular momentum loss due to the stellar wind propagates to the angular
momentum loss of the entire CV system (Verbunt & Zwaan, 1981; van Paradijs,
1986; Hameury et al., 1988; Mestel & Spruit, 1987).

When the orbital period drops to $\sim 3$ hours, the donor star becomes
convective and shrinks to a size below its Roche radius, which disrupts the
magnetic dynamo for sustaining the magnetic braking (Garraffo et al., 2018).
Because the mass loss rate of the donor directly translates to the angular
momentum loss rate, the cessation of magnetic braking leads to a sudden
decline of both these physical quantities. The accretion ceases and the system
evolves as a detached binary, so for a while the orbital period shortens only
due to gravitational radiation, until the donor star re-fills its shrunken
Roche lobe (Spruit & Ritter, 1983). The re-attachment of the donor to the
Roche lobe at around $P_{\rm orb}=2$ hours re-initiates the mass transfer, and
the system re-emerges as an active CV. As a result, there is a paucity of CVs
with periods between 2 hours and 3 hours, a phenomenon known as the period gap
(Rappaport et al., 1982; Verbunt, 1984; Ritter & Kolb, 2003; Knigge, 2006).

The orbital period continues to shorten for these re-ignited CVs via
gravitational radiation and mass loss. The mass transfer eventually causes the
system to go out of thermal equilibrium due to the Kelvin-Helmholtz time scale
growing faster than the mass-loss time scale. In this situation, further mass
loss results in an increase in orbital period (Paczynski & Sienkiewicz, 1981),
and the period distribution of CVs displays a minimum at P${}_{\rm min}\simeq
80$ min (Kolb, 1993; Kolb & Baraffe, 1999; Ritter & Kolb, 2003; Gänsicke et
al., 2009).

In this paper, we present observations and modeling of an enigmatic CV CSS
J160346.08+193540.0 (hereafter CSS1603+19). This object came to our attention
as a periodic infrared variable with a period of P$=81.96$ min – near the
theoretical minimum for CVs – identified during our exploration of the time-
variability data from the Wide-field Infrared Survey Explorer (WISE; Wright et
al., 2010) which included the search for eclipsing binaries and other periodic
variables (Hwang & Zakamska, 2020; Petrosky et al., 2021). CSS1603+19 exhibits
a strong mid-infrared excess which cannot be accounted for by a brown dwarf
companion, and a mid-infrared periodic variability of one magnitude. The
source was known previously as a candidate “polar” or AM Herculis star from
its photometric characteristics (Oliveira et al., 2020) and optical
spectroscopy, with the radial velocity variability of emission lines
indicating the same orbital period as the infrared light curve (Thorstensen et
al., 2017). Polars are CVs in which the magnetic field is high enough to
prevent the formation of an accretion disk, to channel the accreting flow into
streams along the field lines and to therefore synchronize the rotation of the
white dwarf with the orbit.

In Section 2, we describe the existing survey and follow-up data of
CSS1603+19. In Section 3 we explain our fitting procedures for photometric,
spectroscopic and light curve data. We present the analysis of these
measurements and a geometric model in Section 4. We conclude in Section 5.
Emission-line wavelengths are given in air. To describe long-term variation
over an extended period (on timescales of months) we use the terminology of
‘high’ vs ‘low’ states, and we use ‘bright’ vs ‘faint’ for the variability
over the 82-minute orbital period.

##  2 Data

In this section, we describe the data we obtained for CSS1603+19. Our data
include both targeted follow-up observations, as well as archival and survey
data for the source.

Figure 1: Left panel: Gaia color-magnitude diagram constructed using sources
within 250 pc from the Sun. The gold-star marker shows the DR3 Gaia average
measurement for CS1603+19. The red dots show the Gaia epoch photometry of
CSS1603+19. Right panel: Raw photometry data, offset vertically for clarity as
listed on the label. Dashed, dotted and solid lines delineate the epochs of
optical MDM, optical APO and NIR GNIRS spectroscopy, correspondingly.

###  2.1 WISE Photometry

CSS1603+19 is well detected by WISE in W1 (3.6 µm) and W2 (4.5 µm) with
signal-to-nose ratios $>28$ in the AllWISE catalog, but is not detected in W3
(12 µm) and W4 (22 µm). CSS1603+19 was covered by both AllWISE and NEOWISE
(Mainzer et al., 2011) eras of the mission. Thus, we were able to obtain
reliable mid-infrared light curves for $\sim 11$ years from 2010 to 2021,
totaling 520 data points for single-epoch photometry.

We used W1 (3.6 µm) and W2 (4.5 µm) band data to construct the photometric
light curve. The data were downloaded from NASA/IPAC Infrared Science Archive.
To ensure the quality of single-epoch exposure, we followed the guideline
provided by Cutri et al. (2012) by adopting the criteria of saa_sep > 0,
moon_masked = 0, qi_factor > 0.9 for
AllWISE111https://wise2.ipac.caltech.edu/docs/release/allwise/expsup/sec3_1.html,
and saa_sep > 0, moon_masked = 0, qi_factor > 0.9, qual_frame > 0 for
NEOWISE222https://wise2.ipac.caltech.edu/docs/release/neowise/expsup/sec3_1a.html.
Redundant exposures contained in AllWISE Multiepoch Photometry Database were
excluded by matching source ID (source_id_mf) in the AllWISE and the NEOWISE
catalogs333https://wise2.ipac.caltech.edu/docs/release/allwise/expsup/sec3_2.html.

WISE photometry is measured in the Vega magnitude system. For consistency in
this paper, we converted the W1 and W2 photometry to the AB-magnitude system
using the following
equations444https://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec4_4h.html:

$\displaystyle W1_{AB}=W1_{Vega}+2.699,$ $\displaystyle
W2_{AB}=W2_{Vega}+3.339.$

###  2.2 ZTF Photometry

We constructed the optical light curve of CSS1603+19 using data from Zwicky
Transient Facility (ZTF; Masci et al., 2018; Dekany et al., 2020) at the
Palomar Mountain in California. We use DR12 from ZTF, which provides
observations of CSS1603+19 between March 2018 and May 2022 for over 4 years
totaling 804 exposures. We used $g$, $r$, and $i$ bands for the light curves.
The data are queried from NASA/IPAC Infrared Science Archive. Quality filters
are already applied to the light curve data produced by the archive.

###  2.3 Gaia Photometry

We supplemented our optical data using Gaia photometry from Gaia Data Release
3 (DR3; Gaia Collaboration et al., 2016, 2022). Gaia performed optical epoch
photometry for CSS1603+19 from 2014 to 2017 totaling 128 data points from the
G, BP, and RP bands. Gaia data are quality-filtered by the
rejected_by_photometry and rejected_by_variability flags. The Gaia photometry
supplements the optical photometry of ZTF observations. The Gaia parallax of
CSS1603+19 is $4.25\pm 0.12$ mas, giving the distance of the source from the
Solar System to be $d=235.3_{-6.4}^{+6.9}$ parsec. The coordinate of
CSS1603+19 from Gaia DR3 is (ra, dec)=(240.94205667792, 19.59439521850) in
degrees at the reference epoch of J2016.0.

In addition, we obtained $G$-band magnitude and $BP-RP$ color of stars within
250 parsec range from the sun to construct the Hertzsprung-Russell (HR)
diagram (Hertzsprung, 1909; Rosenberg, 1910; Russell, 1914). To control the
signal-to-noise (S/N) ratio of data, we adopted the following criteria for our
query:

parallax_over_error > 10

phot_g_mean_flux_over_error > 50

phot_bp_mean_flux_over_error > 20

phot_rp_mean_flux_over_error > 20

phot_bp_rp_excess_factor < 2

ruwe < 1.2

These criteria yield approximately two million Gaia sources for the color-
magnitude diagram.

###  2.4 GNIRS Near-Infrared Spectroscopy

We obtained the near-infrared (NIR) spectroscopy of CSS1603+19 using the
Gemini Near-Infrared Spectrograph (GNIRS; Elias et al., 2006b; Elias et al.,
2006a) instrument at the Gemini North Observatory on Mauna Kea, Hawaii
(program ID GN-2019A-Q-322, PI: Hwang). We used the cross-dispersed mode with
a 0.3-arcsecond slit, the 32 l mm-1 grating, and the short blue camera (15$"$
pixel-1) for the observation. The setting ensured the S/N ratio is above 3,
and the spectral resolving power is above 1600 to resolve emission lines and
their radial velocity variations. On the nights of 2019-02-13, 2019-03-09, and
2019-03-17, we obtained four, nineteen, and fourteen 300-second exposures
respectively to cover the full orbital period of CSS1603+19. Our NIR spectrum
covers wavelength range between 8500 Å and 26000 Å, giving full coverage over
the $J$, $H$, and $K$ bands.

We reduced the data using PypeIt v.1.8 and v.1.10 (Prochaska et al., 2020).
PypeIt uses A-B image differencing as the first step of sky subtraction, then
improves sky subtraction further in individual 2D A-B difference spectra using
sky apertures. It uses sky lines for wavelength calibration. We flux-calibrate
using bright A0 stars observed right before and right after science exposures.
The telluric absorption correction is computed from the spectrum itself
assuming an underlying polynomial continuum. As a final step after all PypeIt
reductions are done, we correct the wavelength calibration and the timestamps
of the exposures to the heliocentric frame.

###  2.5 MDM Optical Spectroscopy

We obtained optical spectra at the Michigan-Dartmouth-MIT (MDM) Observatory on
Kitt Peak, Arizona (Thorstensen et al., 2017). The observation used the
“Nellie” detector covering wavelength range between 4460Å and 7770Å. Bright O
and B-stars were used to correct for telluric absorption. Comparison lamps and
night sky line [OI]$\lambda$5577Å were used to track night spectrograph
flexure (Thorstensen & Skinner, 2012). A total of 32 600-second exposures were
taken over multiple days in March and June 2017 in a high state. The raw data
were reduced using IRAF routines. One-dimensional spectra were extracted using
optimal-extraction algorithm of Horne (1986).

###  2.6 APO Optical Spectroscopy

We acquired more optical spectroscopy from the Apache Point Observatory (APO)
at Sunspot, New Mexico. The observation was conducted on August 15, 2019 (PI:
Hwang). The weather conditions were clear. We observed the source with Dual
Imaging Spectrograph, where the red camera covers the wavelength from 6000Å to
7000Å. A 1.5-arcsec slit and the R1200 grating were used, with a resolving
power $\sim 1000$. The blue camera suffered from light contamination and did
not provide scientifically useful data. We observed the target for $\sim$1.5
hours (covering about one full orbital period of CSS1603+19), with three
15-minute exposures and five 10-minute exposures bracketed by wavelength
calibration. The data are reduced by the standard IRAF pipeline.

We were hoping to catch the object in a low state and detect the photosphere
of the white dwarf in the APO observations, but caught it in an active (high)
state, and we do not see any white dwarf photospheric features. We detect
H$\alpha$ line and He I emission lines in the APO spectrum, but the following
analysis focuses on the MDM spectra (Sec. 2.5) for their larger wavelength
coverage and a larger number of exposures. The overall behavior (e.g. radial
velocity variations of H$\alpha$) of the APO spectra is consistent with that
seen in the MDM spectra.

##  3 Measurements

This section describes the photometric and spectroscopic measurements for
CSS1603+19. In Section 3.1, we derive the object’s orbital period from the
mid-infrared light curves using WISE data. In Section 3.2 we analyze optical
photometry from ZTF and Gaia. In Sections 3.3 and 3.4 we present the continuum
emission and strong line emission analysis for GNIRS near-infrared data, and
in Section 3.5 we perform similar analysis but in the optical domain. All
phase-resolved plots have the phase aligned with the WISE zero-phase point.

Figure 2: WISE light curve phase-folded with a period of 81.96283 min. Long-
term variations are removed.

###  3.1 Mid-Infrared Photometry

AllWISE and NEOWISE observations span $\sim 11$ years. The observed light
curve reveals variability on a range of timescales. In particular, the light
curve exhibits a long-term variation, with three high phases captured in 11
years where high phases are approximately one magnitude brighter than the low
phase as shown in Figure 1 (right panel).

We search for periodicity after removing the long-term variability.
Specifically, we group WISE photometry into 30-day groups. For each data
point, we subtract the mean magnitude of its corresponding group and add the
mean magnitude of all data points, thus removing the long-term variation from
the WISE photometry. Then, we use the Lomb-Scargle periodogram as implemented
in astropy (Astropy Collaboration et al., 2018) to search for the periodicity
in the mean-calibrated WISE photometric data (Lomb, 1976; Scargle, 1982). A
prominent peak is detected, allowing us to derive the orbital period of
CSS1603+19 of $P_{\rm orb}=0.05691863(25)$ day, in agreement with the period
measurement $P=0.05693(2)$ day from Thorstensen et al. (2017), from the fits
to the H$\alpha$ emission-line velocity. We see no evidence of any change in
the period over the 11-year timescale of WISE observations.

We obtain the phase-folded light curve of CSS1603+19 shown in Figure 2 by
using $P_{\rm orb}=0.05691863$ day. The v-shaped dimming feature has a minimum
at orbital phase $0.701\pm 0.012$. In all subsequent figures we phase-fold
light curves with the same period and phase zero point as the WISE
observations so that the phases are directly comparable. By fitting the shape
of the light curve with parametric functions or splines, we find that the
object spends $45-50\%$ of the orbital period at flux values below the median
flux. The peak-to-peak amplitude of the light curve is $1.1$ mag in W1 and
$1.0$ mag in W2, but given the quality of the data this difference in depth is
not statistically significant. We roughly split the WISE light curve into
‘high’ and ‘low’ states by considering separately the epoch fluxes above and
below the median and we find that the eclipse depth and shape do not depend on
whether the object is in the ‘high’ or ‘low’ state.

###  3.2 Optical Photometry

Figure 1 shows the Gaia color-magnitude diagram. The Gaia sample consists of
2.2 million sources restricted to distances $<250$ pc from the Sun. The gold
star marks the DR3 mean of the measurements for CSS1603+19 at $M_{G}=11.855$
mag and $BP-RP=0.356$ mag. The DR3 measurements are based on Gaia photometric
measurements between Sept 2014 and Apr 2017. This period captures an active
high phase of CSS1603+19 starting June 2016. The active phase causes the epoch
photometry of CSS1603+19 to be brighter and redder than the white dwarf
sequence, as shown by the red dots in Figure 1. The absolute $G$-band
magnitude of CSS1603+19 reaches a minimum at around 12 mag, and when
CSS1603+19 is in this low state no optical variability is detected.

ZTF observatory provides $\sim 3$ years of optical photometry between 2018 and
2022. The observation captured two high states in Jan 2019 and Aug 2019\. The
high states led to $\sim 2.5$ mag variability in the ZTF’s $r$ and $i$ bands
and $\sim 2$ mag variability in the $g$ band as shown in Figure 1 (right
panel). We construct the phase-folded light curve using the ZTF photometric
measurements in the low state. In contrast with mid-infrared measurements, the
ZTF photometric light curve shows no variability at the orbital period.

Figure 3: A typical GNIRS near-infrared exposure of CSS1603+19. The exposure
was obtained on 2019-03-09. Hydrogen Paschen series emission lines are strong.
Brackett series at high order are also visible within the wavelength range. He
II at $\lambda$10094 and He I at $\lambda$10830 are also evident. Broad peaks
of cyclotron emission are visible at long wavelengths (as indicated by the
arrows and shown in the inset). Shaded regions are strongly affected by
telluric absorption that may contain telluric residual after reduction.

###  3.3 Near-Infrared Spectroscopy – Continuum Features

A typical GNIRS spectroscopic exposure is shown in Figure 3. We observe a
time-varying continuum with strong emission lines and additional more subtle
variations.

We fit a low-degree polynomial to the near-infrared continuum, masking the
emission lines and overly noisy data points. We then calculate $JHK$-band
magnitudes from the fits by integrating the model spectral flux density within
each band. We then phase-fold these synthetic continuum light curves using the
WISE-derived period $P_{\rm orb}$ and its phase zero point, so that the phases
in all the phase-folded light curves are directly comparable. The near-
infrared bands show higher variability compared to the mid-infrared range: $J$
and $H$-band light curves both show peak-to-peak amplitude up to $\sim 3$ mag.
The least variable $K$ band has peak-to-peak amplitude $\sim 2.5$ mag.

The continuum color is not constant. By computing the difference between the
synthetic continuum fluxes, we find that the continuum becomes $\sim 0.4$ mag
redder before the flux minima of the $JHK$ bands.

For some near-infrared spectroscopic exposures, we observe spectral wiggles
between $\sim 20,500$Å and $\sim 24,500$Å as shown in Figure 3. These spectral
features are likely cyclotron emission features originating in the strong
magnetic field near the white dwarf surface. Subtracting the continuum to
reveal these features, we find that they vary in intensity throughout the
orbital period (Figure 4), which is common for polars (Schwope & Beuermann,
1990; Schwope et al., 1993). The NIR cyclotron feature is at its maximum
intensity at the orbital phase $\varphi\simeq 0.13\pm 0.15$ and at its minimum
intensity at $\varphi\simeq 0.63\pm 0.15$ – i.e., approximately in phase with
the WISE photometric intensity variation, which has a minimum at
$\varphi=0.701\pm 0.012$.

We recognize that if the donor is a brown dwarf, then its atmosphere may lead
to broad spectral features in the longer wavelength domain in the near
infrared spectroscopy due to the forest of molecular opacity features (Ciardi
et al., 1998). In comparing the continuum-subtracted spectrum of CSS1603+19
with brown dwarf templates from Peterson et al. (2008), we find that the
trough between the broad peak features in CSS1603+19 is inconsistent with the
wavelength centroid of the dominant H2O absorption in brown dwarfs.
Furthermore, neither the overall flux of the features, nor their orbital
evolution are consistent with a brown dwarf origin (Section 4.1). Thus we
conclude that the atmosphere of the donor is unlikely to be responsible for
these features.

Figure 4: Phase-resolved spectroscopy which reveals cyclotron features in near
infrared. Continuum is subtracted. Clear near-infrared broad peaks appear
during the bright state of the WISE light curve centered around 21400Å and
23900Å.

The central wavelength of the $n$-th cyclotron harmonic is

$\lambda_{n}=\frac{2\pi mc^{2}}{eBn},$ (1)

where $m$ and $e$ are electron’s mass and charge and $B$ is the magnetic
field. One feature is centered at $\sim 21,400$Å and another at $\sim
23,900$Å. If these features are due to cyclotron emission, then this gives us
two measurements for two unknowns ($n$ and $B$), allowing us to estimate
$n\simeq 8-10$ and $B\simeq 5-6$ MG. We do not detect any other cyclotron
harmonics at shorter wavelength; our ability to do so may be affected by
telluric corrections. We discuss the cyclotron scenario in more detail in Sec.
4.3.

###  3.4 Near-Infrared Spectroscopy: Emission Lines

Figure 5: Top left: double-Gaussian fit to continuum-subtracted Paschen
$\beta$ emission line in one exposure, with data in black and overall fit in
blue line. The fit for the broad peak is shown with the yellow dot-dashed
curve, and the narrow component with the red dashed curve. Top right: Paschen
$\beta$ orbital evolution and fits on 2019-03-09. From left to right, the
panels are raw data, narrow peak fit, broad peak fit, double peak fit, and
residual significance (maxed at S/N=3), with yellow corresponding to bright
emission and blue corresponding to faint emission. The double peak emission is
most obvious for exposures 8-11 where the broad component presents on the blue
(short wavelength) side of the narrow peak, and for exposure 1-4 where the
broad component is on the side of the narrow peak. Bottom row, from left to
right: phase-resolved spectra of Paschen $\beta$, Paschen $\gamma$, Balmer
$\alpha$, and Balmer $\beta$ respectively, with black corresponding to bright
emission and white corresponding to faint emission. One NIR exposure is
removed due to cosmic ray contamination. Continuum is subtracted. Two emission
components are visible. The narrow component is traced by the red sinusoidal
curve; the yellow curve traces the broad component. The peak wavelengths of
the two peaks obtained from individual fits are marked by points with error
bars. The broad component has a larger velocity semi-amplitude than the narrow
component and precedes the narrow component by approximately 0.2 of the
orbital phase.

All identified strong lines in Figure 3 show similar orbital variations of
their radial velocity profile. After excluding the lines whose fitting quality
is affected by the continuum emission and other emission mechanisms, we model
five emission lines (Paschen $\beta$, $\gamma$, $\delta$, $\epsilon$, and
Brackett $\gamma$) to derive their radial velocities and other kinematic
measures.

We first fit emission lines using a Gaussian profile. The radial velocities
derived from the Gaussian centroid of each peak are noise-weighted for each
exposure. To model the orbital evolution of the emission lines, we then
perform non-linear least squares fit for the noise-weighted radial velocity
using a sinusoidal model in the phase domain:

$\displaystyle V(\varphi\;|\;\gamma,K,\phi)=\gamma+K\sin(2\pi(\varphi-\phi))$
(2)

where $K$ is the velocity semi-amplitude, $\gamma$ is the systematic velocity
and $\phi$ is the overall phase shift. We minimize the number of free
parameters for the model to maximize reliability of the radial velocity
measurement, for example by fixing the orbital period to that derived from
WISE photometry. The model fitting yields semi-amplitude $K=190.1\pm 7.6$ km
s-1, systemic velocity $\gamma=83.3\pm 5.3$ km s-1, and phase of the velocity
node $\phi=0.5046\pm 0.0059$.

A closer look at the emission lines reveals interesting substructure. In
Figure 5, we show the phase-resolved spectra of Paschen emission lines where
we split the orbital period into 15 bins. Aside from the stronger component
traced by the red line, a fainter component is also visible. Accurate
characterization of the faint component needs a high S/N ratio and minimal
effects from nearby emission peaks and continuum noise. Therefore, while we
detect similar double peaked emission from both hydrogen and helium lines, we
only analyze the double-peaked structure of the strongest lines – Paschen
$\beta$ and Paschen $\gamma$ – in the near-infrared domain.

We use a double-Gaussian model to fit the emission lines in all exposures via
least squares optimization, fixing the ratio of the peak amplitudes throughout
the orbit. The time evolution of each component is modeled using the sinusoid
from Equation 2. Velocity dispersion for both components are free to vary in
the time-dependent fit. The double-Gaussian fitting of the emission lines
throughout the orbit demonstrates that the velocity dispersion of the narrow
peak is consistently smaller ($\sim 150$ km s-1) than that of the broad peak
($\sim 350-600$ km s-1). At the spectral resolution $R\gtrsim 1600$, the
emission lines are well spectrally resolved, so these values reflect the
spread of physical velocities of gas responsible for the emission. Figure 5
shows a typical fit for one exposure for the Paschen $\beta$ line and the
overall fit quality for 18 exposures on 2019-03-09. One exposure is removed
due to potential cosmic ray contamination. The results are summarized in Table
1.

The broad component and the narrow component are systematically offset from
one another by about 0.2 in orbital phase. To verify the reliability of this
measurement, we fit the entire ensemble of emission lines in the velocity
phase space using the model in Equation 3. The Equation describes the double
peak combination in given phase bin $\varphi$ is as follows:

$\displaystyle I(v|\varphi)=$
$\displaystyle\sum_{i\in\\{n,b\\}}c_{i}I_{\varphi}\exp\left(-\frac{(v-V(\varphi\;|\;\gamma_{i},K_{i},\phi_{i}))^{2}}{2\sigma_{i}^{2}}\right)$
(3)

where $I_{\varphi},\gamma_{i},K_{i},\phi_{i}$ are parameters need to be
fitted. We use $N_{\rm bin}=15$ bins for constructing the phase-resolved
spectroscopy, therefore we have $N_{\rm bin}+6$ free parameters in total. Here
$i\in\\{n,b\\}$ are the indices distinguishing the narrow and the broad
Gaussian components, and $I_{\varphi}$ is the amplitude of the narrow peak the
given phase bin $\varphi$. Radial velocity
$V(\varphi\;|\;\gamma_{i},K_{i},\phi)$ has the same sinusoidal definition as
in Equation 2 to model the peak position across the entire orbital phase, but
now we are fitting all spectra as an ensemble. $c_{i}$ are constants used for
fixing the amplitude ratio between narrow and broad peaks. We use $c_{n}=1$
and $c_{b}=0.65$ for a reasonable fitting. $\sigma_{i}$ are the widths of the
narrow and broad Gaussian components.

We fix the Gaussian widths by using the median values we obtained using prior
fitting of all spectra independently, which allows us to limit the number of
free parameters without significantly harming the overall fitting quality. On
the contrary, the peak amplitude shows significant variation over the orbital
phase. Fixing the broad amplitude introduces large fitting residuals,
defeating our goal to verify the systematic velocity measure, so we allow it
to vary. The fits for Paschen $\beta$ and Paschen $\gamma$ confirm that the
phase offsets between the narrow and broad components are $0.203(14)$ and
$0.189(22)$ respectively – i.e., close to, but not quite a quarter of the
orbital period, with the phase difference of 0.05 from the quarter of a period
now verified using ensemble fitting at high statistical significance.

The net systematic radial velocities of the broad and narrow components appear
to be different as well. To verify the difference of systematic velocity, we
compare the fitting result with and without requiring $\gamma_{n}=\gamma_{b}$.
For the two Paschen lines, both the Akaike information criterion and the
Bayesian information criterion show significant favor to the model whose
$\gamma_{n}\neq\gamma_{b}$. The broad peaks are systematically blue-shifted by
$212\pm 15$ km s-1 for Paschen $\beta$ and $176\pm 29$ km s-1 for Paschen
$\gamma$ compared to the narrow peaks. There are differences between the
ensemble fitting results and the results from fitting individual exposures
summarized in Table 1. While the ensemble fitting is better for demonstrating
the statistical significance of the phase difference and systematic velocity
difference between the narrow and broad components of the emission lines, our
analysis below is based on the fitting of individual exposures in which fewer
constraints are imposed over the fit.

Table 1: Double-Gaussian emission line fit result. The mask column indicates the orbital phase masked that is not considered during fitting. Component indicates distinguishes the broad (B) and narrow (N) emission component of the emission line. $\sigma$ is the mean of velocity dispersion values obtained in individual fits. Strong Line | Mask | Comp. |  $K$ (km s${}^{-1})$ |  $\gamma$ (km s${}^{-1})$ | $\phi$ |  $\sigma$ (km s-1)
---|---|---|---|---|---|---
Paschen $\beta$ | $0.65<\phi<0.85$ | N | $192.1\pm 7.5$ | $93.0\pm 5.3$ | $0.6172\pm 0.0074$ | $190\pm 69$
B | $289\pm 28$ | $-55\pm 17$ | $0.414\pm 0.011$ | $620\pm 110$
Paschen $\gamma$ | $0.6<\phi<0.8$ | N | $181.9\pm 10.0$ | $87.6\pm 7.7$ | $0.638\pm 0.011$ | $173\pm 74$
B | $251\pm 25$ | $-22\pm 14$ | $0.408\pm 0.011$ | $600\pm 140$
Balmer $\alpha$ | $0.65<\phi<0.75$ | S | $196.3\pm 6.1$ | $9.4\pm 4.4$ | $0.6433\pm 0.0053$ | $154\pm 36$
W | $318\pm 21$ | $47\pm 13$ | $0.4438\pm 0.0080$ | $360\pm 73$
Balmer $\beta$ | $0.85<\phi<0.9$ | S | $189.7\pm 6.4$ | $37.9\pm 5.3$ | $0.6640\pm 0.0069$ | $169\pm 38$
W | $281\pm 17$ | $15\pm 11$ | $0.4655\pm 0.0070$ | $365\pm 81$

###  3.5 Optical Spectroscopy

The optical spectrum shows clear continuum emission and strong emission lines
as exhibited in Figure 6. Using the same approach used for near-infrared
analysis, we fit a low-degree polynomial to approximate the continuum emission
while masking the emission lines. The continuum emission is more stable over
the optical range compared to its near-infrared counterpart. The continuum
flux remained approximately constant over the orbital period. No significant
pattern is recognized in the synthetic $V$-band and $R$-band light curves.

Figure 6: Mean optical spectra: black for the March 2017 MDM observation, gray
for the June 2017 MDM observation and blue for the August 2019 APO
observation. Balmer series, He I, and He II $\lambda$4686Å are clear, but the
He II line becomes undetectable during the June 2017 observation. The
continuum emission is also fainter in June 2017.

In addition to the strong Balmer series of hydrogen, He I lines and He II
$\lambda$4686Å are evident in the optical range during the March 2017
observation. The He II line becomes undetectable in June 2017 when CSS1603+19
enters a low phase. Even in the faintest optical spectrum we have obtained,
the white dwarf photosphere is not discernible, unlike the case of some CVs
with extremely low accretion rates (Gänsicke et al., 2009).

The two-component emission-line structure is also observed in the optical
spectrum. We first fit for the overall orbital evolution of the radial
velocity without deblending the two components using the model described by
Equation 2. We include Balmer $\alpha,\beta$, He I $\lambda$5016Å, and He I
$\lambda$6678Å. We find semi-amplitude $K=197.0\pm 5.2$ km s-1, systemic shift
$\gamma=35.4\pm 3.4$ km s-1, and phase of the velocity node $\phi=0.4549\pm
0.0041$. The semi-amplitude is slightly larger than the measurement using
near-infrared data, but the difference is within the uncertainty of the
measurement.

We then deblend the broad and the narrow components of emission lines using
the same method as in Section 3.4, fitting Balmer $\alpha$ and Balmer $\beta$.
The fitting result is summarized in Table 1 and visualized in Figure 5.
Similar to the near-infrared case, the broad emission peak precedes the narrow
emission peak by 0.2 orbital phase up to a high accuracy. In contrast with
Paschen $\beta$ and Paschen $\gamma$, the broad components of Balmer lines do
not have a blue systematic shift. By fitting the ensemble of spectra using
Equation 3, we do not find any evidence for a systematic offset between broad
and narrow components of Balmer $\alpha$. The Balmer $\beta$ line has its
broad component systematically redshifted by $37\pm 9.0$ km s-1 compared to
the narrow component. This is in contrast to the systematic blueshift we see
for the Paschen sequence.

##  4 Analysis and Discussion

This section focuses on the physical interpretation and analysis of
measurements from Section 3.

###  4.1 Spectral Energy Distribution (SED) Analysis

Using optical and infrared data from WISE and ZTF and our own spectroscopic
measurements from MDM and Gemini, we construct the spectral energy
distribution (SED) of CSS1603+19 shown in Figure 7. We select available
observations that fall into a relatively bright orbital phase
($0.15<\varphi<0.25$) and into a relatively faint orbital phase
($0.65<\varphi<0.75$). We observed significant variability in the infrared. On
the contrary, the variability in the optical is uncorrelated with the orbital
phase, and this is also apparent from Figure 7. The SED shows a clear infrared
excess and a local minimum at around $\lambda\simeq 8000$Å.

Figure 7: SED of CSS1603+19 in the bright phase (orbital phase bin
$0.15<\varphi<0.25$, top) and in the faint phase ($0.65<\varphi<0.75$,
bottom). The colored dots are observations labeled as by legends. The blue
curve shows a 13,000 K white dwarf model with surface gravity $\log g=7.75$
from Koester (2010) and Tremblay & Bergeron (2009) as our white dwarf
estimate. The orange and gold curve showed the 0.1 R⊙ brown dwarf emission
model at 1600K and 1300K respectively using templates from Marley et al.
(2021).

We first attempt to model the SED of CSS1603+19 with a sum of spectra of a
white dwarf and a brown dwarf. We select a $13000K$, $\log g=7.75$ white dwarf
template from Koester (2010) and Tremblay & Bergeron (2009) to model the white
dwarf emission. In some cases, the brown dwarf atmospheric absorption bands
may lead to spectral features between 2 µm and 2.4 µm (Ciardi et al., 1998)
qualitatively similar to those seen in our data. Thus, we use brown dwarf
atmosphere templates from Marley et al. (2021) spanning a temperature range of
800K to 2400K with surface gravity $\log g=5$ to check whether our spectral
wiggles match the brown dwarf absorption features. The brown dwarf H2O
absorption leads to a broad peak between 2.1 µm and 2.3 µm whose exact
centroid is dependent on the brown dwarf temperature. To match the peak at
2.14 µm observed in our target, the brown dwarf temperature needs to be
between the 1300 and 1600 K range.

Upon determining the temperature, we normalize the spectrum using the brown
dwarf radius $0.1R_{\odot}$, leading to the results shown in figure 7. By
comparing the model spectral peak and the observed spectral peak, we find that
the flux of the peak in the 1600 K 0.1 $R_{\odot}$ model is too faint by a
factor of 10 compared to the observed value, and the flux in the 1300K model
is too weak by a factor of 30. Scaling up the brown dwarf to match the model
flux within the peak to the observed one yields a minimum brown dwarf radius
$\sim$0.35 $R_{\odot}$, which is nonphysical. Apart from the peak at 2.14 µm,
the brown dwarf model also predicts two other peaks near 1.65 µm and 1.3 µm,
which are undetected in our NIR spectroscopy, while the observed peak at 2.39
µm peak in our spectroscopy remains unaccounted for by the brown dwarf
template.

The overall luminosity of our model brown dwarfs is too faint by about two
orders of magnitude compared to the observed infrared continuum in the bright
phase. In CVs, the donor star is significantly larger than the white dwarf.
Thus, the large periodic variability of the 2.14 µm feature and the NIR
continuum are also unexplained in a model where the NIR flux is dominated by
the brown dwarf, as there is no occulting body with the size comparable to the
size of the brown dwarf.

In some CVs, there is a circumbinary disk that can produce infrared excess.
The size of the circumbinary disk may help explain the strong infrared
component in the spectral energy distribution which cannot be accounted for by
the donor star (e.g. Dubus et al. 2004). However, the geometry of the
circumbinary disk is not compatible with the eclipsing-like infrared
lightcurve and rapid variability. The infrared emission from a circumbinary
disk can only be minimally blocked by the binary, and the geometry prevents
the observed factor of two variability. Thus, the model cannot explain the
large 1-mag infrared variability in CSS1603+19. Several CVs suspected of
having circumbinary disks were later shown to have their infrared excess
dominated by the donor or by cyclotron emission (Harrison et al., 2013; Dubus
et al., 2004).

###  4.2 Binary model

Figure 8: The geometry of CSS1603+19. The observer is fixed below the figure
with the line of sight indicated by the arrow, while the system rotates
counterclockwise. Gold and red dots represent the line-emitting regions
producing the broad and narrow components, respectively. The same color scheme
is used in Figure 5 to illustrate the radial velocity behavior of these lines.
The green spot marks the location where the stream’s geometry becomes shaped
by the white dwarf magnetic field. $\varphi^{\prime}=0.6874\pm 0.0089$ is the
orbital phase where the broad component appears with the largest redshift.
$\varphi^{\prime\prime}=0.8915\pm 0.0071$ is the orbital phase where the
narrow component appears with the largest redshift. The proposed geometry of
the curved accretion stream results in
$\varphi^{\prime\prime}-\varphi^{\prime}=0.204\pm 0.016$ (instead of the 0.25
if the stream were straight). While the white dwarf is also moving relative to
the center of mass of the system, its mass is likely much larger than that of
the donor, and therefore we illustrate it as nearly stationary. The orbital
inclination is not incorporated in the figure. Table 2: Stellar Parameters
System Parameters
---
$P=81.96282\;\pm\;0.00036\textrm{\;min}$
$d=235.3_{-6.4}^{+6.9}\textrm{\;pc}$
$i\lesssim 74.6\pm 0.4^{\circ}$
$q\lesssim 0.263\pm 0.046$
$a\gtrsim 0.425\pm 0.024R_{\odot}$
Stellar Mass Radius
$M_{1}\gtrsim 0.245\pm 0.032M_{\odot}$
$R_{1}\lesssim 0.01241\pm 0.00023R_{\odot}$
$M_{2}=0.0644\pm 0.0074M_{\odot}$
$R_{2}=0.113\pm 0.013R_{\odot}$
Cyclotron Region
$B=5-6MG$

In this section, we constrain the stellar parameters of our binary system and
model a range of observables we discussed above. The first striking feature
that gives us a clue as to the geometry of the system is the radial velocity
variations of the two peaks of the emission lines which are out of phase and
have different amplitudes. Furthermore, the semi-amplitude of the broad
component $K_{b}\simeq 300$ km s-1 is high enough that the locus of the
possible solutions in the space of white dwarf ($M_{1}$) and donor companion
($M_{2}$) masses suggests unrealistically high masses for the companion ($\ga
0.5M_{\odot}$) which we would have seen in our data. This observation is
reminiscent of that by Mason et al. (2019) for eclipsing polar CRTS J0350+3232
with an orbital period of 2.4 hours. We therefore adapt the geometric model of
Mason et al. (2019) for our system, albeit with some modifications as
described below.

Specifically, following Mason et al. (2019) we hypothesize that the broad
component of the emission line originates from the accretion spot close to the
surface of the white dwarf. Its large radial velocity amplitude is due to a
net flow of matter within the accretion stream toward the white dwarf. Figure
8 illustrates that $\varphi^{\prime}$ is the orbital phase where the radial
velocity of the broad component (gold spot in the figure) is at its maximal
redshift. The narrow component (red dot in Figure 8) is placed nearer to the
donor, although following the investigation of the resulting $M_{1}-M_{2}$
solutions we are inclined to place it at the Roche point rather than at the
exact orbital location of the donor.

If the stream connecting the donor and the white dwarfs was straight, then
just by associating one end of it with the Roche point and the other end of it
with the accretion spot we would expect an 0.25 phase offset between the two
components. In practice, due to the Coriolis force, the stream bends in the
same direction as the orbital motion, and the resulting trajectory was
computed by Lubow & Shu (1975). Associating one end of such stream with the
Roche point and the other end with the accretion spot, one would expect an
$>0.25$ phase offset between two components assuming the absence of an
accretion disk. Instead the phase offset is
$\varphi^{\prime\prime}-\varphi^{\prime}=0.204\pm 0.016$. This suggests that
the accreting stream lags the donor star by $\sim 18$ degrees (the same lag is
seen in the polars described by Mason et al. 2019 and Tovmassian et al. 1997).
The bending of the stream required by these observations is possible for a
polar where the accretion geometry close to the white dwarf is dominated by
its magnetic field. Following the AM UMa model depicted in Bonnet-Bidaud et
al. (1996), we show in Figure 8 that as the stream propagates from the donor
to the white dwarf, it first curves in the direction predicted by Lubow & Shu
(1975) and then twists due to the effects of the magnetic field to account for
the observed phase differences in the emission lines. Finally, the $\sim 100$
km s-1 systemic blueshift of the broad component relative to the narrow
component seen in some of the lines (Pa $\alpha$, Pa $\beta$ as listed in
Table 1) may be due to partial obscuration of the accretion spot by the stream
at phase $\varphi^{\prime}$, when the emitting component is at its highest
redshift.

The possible cyclotron features in the near-infrared spectroscopy indicate the
presence of strong magnetic fields from the system, likely associated with the
white dwarf. The strong magnetic field near the white dwarf may interrupt the
accretion disk, causing the accretion stream to flow along the magnetic field
lines. The presence of stable weak broad components to the emission lines that
show regular orbital variations implies a stable accretion flow near the white
dwarf. This indicates the white dwarf rotates synchronously with respect to
the system, revealing the polar (AM Her star) nature of the system in contrast
to an intermediate polar (DQ Her Type) in which the white dwarf rotates
asynchronously (Cropper, 1990). Abundant emission lines and the strong He II
at $4686$Å also mark the system as an AM Her type CV, although the system does
not quite meet the HeII$4686$Å/H$\beta>0.4$ criterion considered sufficient
for classifying a system as a magnetic CV (Silber, 1992; Mason et al., 2019).

Our first dynamical constraint comes from the size of the Roche lobe. The
$81.96$-min orbital period of CSS1603+19 suggests that the system is below the
CV period gap. As accretion signatures are clearly visible, the donor must be
overflowing its Roche lobe whose radius is completely determined by the semi-
major axis and the mass ratio of the system (Eggleton, 1983). Fully convective
stars – such as the donor below the period gap – are well described by $n=1.5$
polytropes. Thus, the donor radius can be constrained using the donor mass-
radius relation and the Roche-lobe radius approximation:

$\displaystyle\frac{R_{2}}{a}=\frac{0.5126q^{0.7388}}{0.6710q^{0.7349}+\ln(1+q^{0.3983})},$
(4)

in which $R_{2}$ is the donor radius, $a$ is the the semi-major axis, and
$q=M_{2}/M_{1}$ is the mass ratio (Sirotkin & Kim, 2009; Knigge et al., 2011).
If we instead assume that the system is a period bouncer, then the period
bouncer mass-radius relation from Knigge et al. (2011) yields a donor mass of
0.1 $M_{\odot}$, somewhat outside of the range of applicability of the mass-
radius relation for bouncers ($M<0.069M_{\odot}$), and we therefore do not
obtain a self-consistent solution. We thus conclude that our target has likely
not yet entered the period-bouncer stage of its evolution; however, our
derived mass is close to the bouncer range and therefore this conclusion is
somewhat uncertain.

We use the power-law mass-radius relation from Knigge et al. (2011) for the
donors of short-period CVs:

$\displaystyle\frac{R_{2}}{R_{\odot}}=0.225\pm
0.008\left(\frac{M_{2}}{M_{\textrm{conv}}}\right)^{0.61\pm 0.01},$ (5)

where $M_{\textrm{conv}}=0.20\pm 0.02M_{\odot}$ is the donor mass at the
period gap (where the donor is fully convective). In contrast to the ordinary
brown dwarf mass-radius relation, the relation is constrained by the period-
density relation of Roche-lobe filling stars on both edges of the period gap
(Knigge, 2006). Finally, for Keplerian orbits

$a=\left(\frac{P_{\rm orb}^{2}G(M_{1}+M_{2})}{4\pi^{2}}\right)^{1/3}.$ (6)

Plugging Equations 5 and 6 into Equation 4, we obtain one equation for two
variables, $M_{1}$ and $M_{2}$. The locus of the resulting solutions in the
$M_{1}-M_{2}$ space (solid line in middle panel of Figure 9) yields
$M_{2}=0.0644\pm 0.0074M_{\odot}$ which is nearly independent of $M_{1}$. The
quoted uncertainty incorporates the stated uncertainties in Equation 5, as
well as marginalizing the result over all possible values of $M_{1}$. This
donor mass is on the boundary between low-mass CVs and period bouncers by the
criteria from Knigge et al. (2011) as expected for the system approximately at
the theoretical CV period minimum.

Our second constraint on the dynamics of the system comes from the radial
velocity measurements of the emission lines. The broad component likely has a
net streaming velocity and therefore cannot be used for determining the masses
of the objects without introducing additional parameters for the motion of the
stream. Placing the narrow component at the center of mass of the donor
results in a solution with an overly small and statistically unlikely
inclination or an unphysically heavy donor. Therefore, we place the emission
region producing the narrow component at the Roche point, following Bonnet-
Bidaud et al. (1996), and we assume that the gas is overflowing the Roche
point with a negligible net velocity, so the observed velocities are entirely
due to the orbital motion at this position. This placement of the narrow
component would result in its total or partial blocking by the donor (i.e.,
dimming) during the phases when the donor is closest to the observer
(analogous to HS Cam as deduced by Tovmassian et al. 1997), which is in
agreement with what we see in Figure 5. We then use the modified binary mass
function to construct a relation between the white dwarf mass $M_{1}$, donor
mass $M_{2}$ and system orbital inclination $i$ (defined relative to the plane
of the sky):

$\displaystyle\frac{M_{1}}{(M_{1}+M_{2})^{\frac{2}{3}}}=(M_{1}+M_{2})^{\frac{1}{3}}\frac{R_{2}}{a}\left(q\right)+\left(\frac{K_{RL}^{3}P_{\rm
orb}}{2\pi G\sin^{3}i}\right)^{\frac{1}{3}}$ (7)

where $K_{RL}$ is the amplitude of the radial velocity of the narrow component
which we assume is located at the Roche point and displays purely orbital
motion, and $\frac{R_{2}}{a}\left(q\right)$ is the Roche lobe radius as a
fraction of the semi-major axis from Equation 4.

The constraints are visualized in Figure 9. As we do not detect a flat-
bottomed photometric dip in the near infrared lightcurve, nor any kind of
periodic variability in the optical lightcurve, we classify the system as non-
eclipsing. Given the size of the donor, this results in an inclination upper
limit of $i\sim 75^{\circ}$, constraining the white dwarf mass to be
$M_{1}\gtrsim 0.262M_{\odot}$. Unfortunately, the orbital evolution of the the
broad line and of the IR photometric light curve are likely strongly shaped by
the magnetic field geometry of the white dwarf, which can have any orientation
relative to the orbit, and therefore we cannot use these observations to place
any other constraints on the orbital inclination. The best estimate of the
binary model parameters are listed in Table 2.

Figure 9: Left: the dynamical constraints on the orbital parameters. The Roche
constraint (that the donor must be at the Roche limit for the system),
combined with a donor mass-radius relationship yields the nearly vertical
black line, i.e., a nearly definite mass for the donor. Each color represents
the $M_{1}-M_{2}$ constraint from the orbital motion of the narrow component
at a fixed inclination. Right: the relationship between the white dwarf mass
and system inclination. Given the mass of the donor, the mass of the white
dwarf and the inclination are constrained by the kinematic solution to lie on
the solid black line. The dot-dashed line marks the upper limit on inclination
to avoid a full eclipse of the white dwarf. The allowed parameter space is of
the system is shaded in grey. The two curved lines represent the solutions of
Equation 7 at the extreme ends of uncertainties in $M_{2}$, $R_{2}$ and
$K_{RL}$.

###  4.3 The origin of infrared emission and variability

The infrared excess of CSS1603+19 at 1-5 µm is too bright to be explained by
the donor alone. Its short-term variability at the orbital period is
surprisingly stable over the 11-year observations from WISE (Fig. 2), with a
peak-to-peak amplitude of $\sim 1$ mag. While there is no significant
difference in the amplitude between W1 and W2, we find a strong wavelength
dependence of the infrared variability in the wavelength range of GNIRS
spectra between 1-2 µm, with the shorter wavelengths having larger amplitude.
In addition to the smooth continuum in GNIRS spectra, there are wiggle-like
spectral features peaking around 2.14 and 2.39 µm (Fig. 4). Such short-term
periodic variability is only present in infrared but not in optical.
CSS1603+19 is among the reddest polars in terms of the cataloged G-W1 color
(Ritter & Kolb, 2003). Its 1-mag peak-to-peak amplitude in WISE is also among
the largest ones in the polars where their WISE light curves are studied in
Harrison & Campbell (2015).

We have ruled out the donor as the dominant contributor to the infrared
emission. In Sec. 4.2, we show that the donor has a radius of $\sim
0.1R_{\odot}$ according to the observed mass-radius relation of donors (Knigge
et al., 2011). If the continuum in GINRS were due to thermal radiation, its
peak at 1.1 µm in Fig. 7 suggests a thermal radiation with a temperature of
$\sim 2600$ K, but the observed infrared luminosity is one order of magnitude
larger than a $0.1R_{\odot}$ thermal source. Therefore, the infrared excess is
not dominated by the donor. This argument also rules out infrared emission due
to the accretion-irradiated donor (Tovmassian et al., 1997): although the
time-variability of the infrared emission is consistent with the self-eclipse
of the donor’s accretor-facing side by the donor itself, the size of required
emitting region would still be well above $0.1R_{\odot}$ and would be
incompatible with our dynamical model. The strong magnetic fields of polars
prohibit the formation of an accretion disk surrounding the white dwarf, and
thus the infrared emission is not from a disk.

Infrared excess and infrared variability at orbital periods in polars are
often associated with cyclotron emission, a byproduct of the accretion shock
at the polar caps of the magnetic white dwarf (e.g. Cropper et al., 1990;
Ferrario et al., 1996; Debes et al., 2006; Campbell et al., 2008a, b; Campbell
et al., 2008c; Harrison & Campbell, 2015). Cyclotron emission can have
multiple harmonic frequencies that span a wide wavelength range, and the
wavelengths of the harmonics depend on the strength of magnetic fields.
Although the exact wavelengths of the harmonics are uncertain due to the broad
spectral profiles of the wiggles, the magnetic field strength needs to be
5-5.5 MG so that $n=8,9$ harmonics match to the observed wavelengths of 2.14
and 2.39 µm, (Sec. 3.3). With a magnetic field strength of $5$ MG, its $n=5,6$
cyclotron harmonics would peak at 3.57 and 4.28 µm, which may explain the
variability in W1 (3.4 µm) and W2 (4.5 µm).

The origin of the infrared variability at its orbital period remains a
mystery. The variability is not due tidal deformation (ellipsoidal variation)
of the donor star because ellipsoidal variation would have a variability
period at half of the orbital period with an amplitude of $\sim 0.3$ mag at
most (Russell, 1945; McClintock et al., 1983), and the donor flux is
insufficient to explain the observed excess for reasonable donor models anyway
(Figure 7). If the infrared emission is dominated by cyclotron near the
surface of the white dwarf, the emitting region may be blocked by the white
dwarf itself due to the white dwarf’s spin (so-called ‘self-eclipse’). The
self-eclipse scenario can naturally explain the stable infrared periodic
variability over 11 years of the WISE data. Since the polar caps are much
smaller than the size of the white dwarf, the v-shaped light curves of WISE
(Fig. 2) require partial self-eclipses of the polar cap; otherwise, a long
duration of full self-eclipses would result in a flat feature in the faint
phase of the light curve. The polar cap eclipsed by the donor star is unlikely
because that would also induce strong variability in optical at the orbital
period, inconsistent with our results.

As seen in Fig. 4, the 2.14 and 2.39 µm spectral features are at their flux
minimum right at the time when the broad emission-line component is at its
maximum redshift (phase $\varphi^{\prime}=0.69$ in Fig. 8). Therefore, it is
the back side of the white dwarf – the one not facing the stream and the donor
– which is producing cyclotron emission. The comparison with emission lines
provides key information about the emitting location of the 2.14 and 2.39 µm
spectral features.

Although the self-eclipsed cyclotron scenario may be able to explain the
infrared excess and variability, there are other observations that may
challenge this picture. In particular, most of the observed polars have
magnetic fields $>10$ MG, and higher-order ($n\geq 6$) harmonics are usually
weak. Hence, the $n=8,9$ harmonics with the magnetic field strength of 5-5.5
MG may not have cyclotron emission as strong as required by the data. It is
possible that the observed infrared variability is a more complex interplay
involving other structures (e.g. eclipse between accretion streams, hot spots,
and/or polar caps). Future X-ray observations may be needed to improve the
geometric constraints on the polar caps (e.g. Heise et al. 1985) to decipher
the mysterious infrared variability.

##  5 Conclusions

In this paper we have presented a follow-up study of a cataclysmic variable
CSS1603+19, which we categorize as a low-mass AM Herculis star. The white
dwarf accretor is likely magnetized, leading to cyclotron emission which
dominates emission at mid-infrared wavelengths. We tentatively measure
cyclotron harmonics in its near-infrared spectrum, allowing us to estimate the
magnetic field strength at $\sim 5MG$.

Using NIR and optical spectroscopy, we detect double-peaked emission lines
whose velocities and amplitudes vary throughout the orbital period. We build a
geometric model of the system (Figure 8) in which the broad component
originates from the accretion stream near the white dwarf and the narrow
component originates from the Roche lobe near the Lagrange point. This model
well explains the orbital variability of the lines, and if we further allow
for the stream to be bent by the white dwarf’s magnetic field away from the
direction implied by the orbital dynamics of the gas, it explains the phase
offset between the components as well.

Requiring that the donor must fill its Roche lobe and that it obeys the
standard mass-radius relationship for the CV donors, we constrain the mass of
the donor to be $M_{2}=0.0644\pm 0.0074M_{\odot}$, with an accuracy of about
10% (this uncertainty is dominated by the uncertainty in the mass-radius
relationship we use). The kinematics of the emission lines and the lack of the
optical eclipse of the accreting white dwarf then allows us to place
constraints on the white dwarf mass ($M_{1}>0.24M_{\odot}$) and the
inclination ($i\lesssim 75^{\rm o}$). The system is likely caught right before
its period bounce.

## Acknowledgments

Y.L., H.-C.H. and N.L.Z. acknowledge support of NASA ADAP grant 80NSSC19K0581
and support from JHU via the seed grant from the Institute for Data Intensive
Engineering and Science. H.-C.H. acknowledges the support of the Infosys
Membership at the Institute for Advanced Study. N.L.Z. was supported at the
IAS by the J. Robert Oppenheimer Visiting Professorship and the Bershadsky
Fund.

This work is based on observations obtained at the international Gemini
Observatory, a program of NSF’s NOIRLab, which is managed by the Association
of Universities for Research in Astronomy (AURA) under a cooperative agreement
with the National Science Foundation on behalf of the Gemini Observatory
partnership: the National Science Foundation (United States), National
Research Council (Canada), Agencia Nacional de Investigación y Desarrollo
(Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina),
Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and
Korea Astronomy and Space Science Institute (Republic of Korea).

This publication makes use of data products from the Wide-field Infrared
Survey Explorer, which is a joint project of the University of California, Los
Angeles, and the Jet Propulsion Laboratory/California Institute of Technology,
funded by the National Aeronautics and Space Administration.

This work mades use of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

This work is based on observations obtained with the Samuel Oschin Telescope
48-inch and the 60-inch Telescope at the Palomar Observatory as part of the
Zwicky Transient Facility project. ZTF is supported by the National Science
Foundation under Grants No. AST-1440341 and AST-2034437 and a collaboration
including current partners Caltech, IPAC, the Weizmann Institute for Science,
the Oskar Klein Center at Stockholm University, the University of Maryland,
Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium
of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin,
Lawrence Livermore National Laboratories, IN2P3, University of Warwick, Ruhr
University Bochum, Northwestern University and former partners the University
of Washington, Los Alamos National Laboratories, and Lawrence Berkeley
National Laboratories. Operations are conducted by COO, IPAC, and UW.

This work is based on observations obtained at the MDM Observatory, operated
by Dartmouth College, Columbia University, Ohio State University, Ohio
University, and the University of Michigan.

This work is based on observations obtained with the Apache Point Observatory
3.5-meter telescope, which is owned and operated by the Astrophysical Research
Consortium.

This research has made use of the NASA/IPAC Infrared Science Archive, which is
funded by the National Aeronautics and Space Administration and operated by
the California Institute of Technology.

## References

  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Bonnet-Bidaud et al. (1996) Bonnet-Bidaud J. M., Mouchet M., Somova T. A., Somov N. N., 1996, A&A, 306, 199
  * Campbell et al. (2008a) Campbell R. K., Harrison T. E., Schwope A. D., Howell S. B., 2008a, ApJ, 672, 531
  * Campbell et al. (2008b) Campbell R. K., Harrison T. E., Mason E., Howell S., Schwope A. D., 2008b, ApJ, 678, 1304
  * Campbell et al. (2008c) Campbell R. K., Harrison T. E., Kafka S., 2008c, ApJ, 683, 409
  * Ciardi et al. (1998) Ciardi D. R., Howell S. B., Hauschildt P. H., Allard F., 1998, ApJ, 504, 450
  * Cropper (1990) Cropper M., 1990, Space Sci. Rev., 54, 195
  * Cropper et al. (1990) Cropper M., et al., 1990, MNRAS, 245, 760
  * Cutri et al. (2012) Cutri R. M., et al., 2012, Explanatory Supplement to the WISE All-Sky Data Release Products, Explanatory Supplement to the WISE All-Sky Data Release Products
  * Debes et al. (2006) Debes J. H., López-Morales M., Bonanos A. Z., Weinberger A. J., 2006, ApJ, 647, L147
  * Dekany et al. (2020) Dekany R., et al., 2020, PASP, 132, 038001
  * Dubus et al. (2004) Dubus G., Campbell R., Kern B., Taam R. E., Spruit H. C., 2004, MNRAS, 349, 869
  * Eggleton (1983) Eggleton P. P., 1983, ApJ, 268, 368
  * Elias et al. (2006a) Elias J. H., Rodgers B., Joyce R. R., Lazo M., Doppmann G., Winge C., Rodríguez-Ardila A., 2006a, in McLean I. S., Iye M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 6269, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. p. 626914, doi:10.1117/12.671765
  * Elias et al. (2006b) Elias J. H., Joyce R. R., Liang M., Muller G. P., Hileman E. A., George J. R., 2006b, in McLean I. S., Iye M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 6269, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. p. 62694C, doi:10.1117/12.671817
  * Ferrario et al. (1996) Ferrario L., Bailey J., Wickramasinghe D., 1996, MNRAS, 282, 218
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2022) Gaia Collaboration et al., 2022, arXiv e-prints, p. arXiv:2208.00211
  * Gänsicke et al. (2009) Gänsicke B. T., et al., 2009, MNRAS, 397, 2170
  * Garraffo et al. (2018) Garraffo C., Drake J. J., Alvarado-Gomez J. D., Moschou S. P., Cohen O., 2018, The Astrophysical Journal, 868, 60
  * Hameury et al. (1988) Hameury J. M., King A. R., Lasota J. P., Ritter H., 1988, MNRAS, 231, 535
  * Harrison & Campbell (2015) Harrison T. E., Campbell R. K., 2015, ApJS, 219, 32
  * Harrison et al. (2013) Harrison T. E., Hamilton R. T., Tappert C., Hoffman D. I., Campbell R. K., 2013, AJ, 145, 19
  * Heise et al. (1985) Heise J., Brinkman A. C., Gronenschild E., Watson M., King A. R., Stella L., Kieboom K., 1985, A&A, 148, L14
  * Hertzsprung (1909) Hertzsprung E., 1909, Astronomische Nachrichten, 179, 373
  * Horne (1986) Horne K., 1986, PASP, 98, 609
  * Hwang & Zakamska (2020) Hwang H.-C., Zakamska N. L., 2020, MNRAS, 493, 2271
  * Kalomeni et al. (2016) Kalomeni B., Nelson L., Rappaport S., Molnar M., Quintin J., Yakut K., 2016, ApJ, 833, 83
  * Knigge (2006) Knigge C., 2006, MNRAS, 373, 484
  * Knigge et al. (2011) Knigge C., Baraffe I., Patterson J., 2011, ApJS, 194, 28
  * Koester (2010) Koester D., 2010, Mem. Soc. Astron. Italiana, 81, 921
  * Kolb (1993) Kolb U., 1993, A&A, 271, 149
  * Kolb & Baraffe (1999) Kolb U., Baraffe I., 1999, Monthly Notices of the Royal Astronomical Society, 309, 1034
  * Lomb (1976) Lomb N. R., 1976, Ap&SS, 39, 447
  * Lubow & Shu (1975) Lubow S. H., Shu F. H., 1975, ApJ, 198, 383
  * Mainzer et al. (2011) Mainzer A., et al., 2011, ApJ, 731, 53
  * Marley et al. (2021) Marley M. S., et al., 2021, ApJ, 920, 85
  * Masci et al. (2018) Masci F. J., et al., 2018, Publications of the Astronomical Society of the Pacific, 131, 018003
  * Mason et al. (2019) Mason P. A., Wells N. K., Motsoaledi M., Szkody P., Gonzalez E., 2019, MNRAS, 488, 2881
  * McClintock et al. (1983) McClintock J. E., Petro L. D., Remillard R. A., Ricker G. R., 1983, ApJ, 266, L27
  * Mestel & Spruit (1987) Mestel L., Spruit H. C., 1987, MNRAS, 226, 57
  * Oliveira et al. (2020) Oliveira A. S., Rodrigues C. V., Martins M., Palhares M. S., Silva K. M. G., Lima I. J., Jablonski F. J., 2020, AJ, 159, 114
  * Paczynski & Sienkiewicz (1981) Paczynski B., Sienkiewicz R., 1981, ApJ, 248, L27
  * Patterson (1984) Patterson J., 1984, ApJS, 54, 443
  * Peterson et al. (2008) Peterson D. E., et al., 2008, ApJ, 685, 313
  * Petrosky et al. (2021) Petrosky E., Hwang H.-C., Zakamska N. L., Chandra V., Hill M. J., 2021, MNRAS, 503, 3975
  * Prochaska et al. (2020) Prochaska J. X., et al., 2020, Journal of Open Source Software, 5, 2308
  * Rappaport et al. (1982) Rappaport S., Joss P. C., Webbink R. F., 1982, ApJ, 254, 616
  * Rappaport et al. (1983) Rappaport S., Verbunt F., Joss P. C., 1983, ApJ, 275, 713
  * Ritter & Kolb (2003) Ritter H., Kolb U., 2003, A&A, 404, 301
  * Rosenberg (1910) Rosenberg H., 1910, Astronomische Nachrichten, 186, 71
  * Russell (1914) Russell H. N., 1914, Popular Astronomy, 22, 275
  * Russell (1945) Russell H. N., 1945, ApJ, 102, 1
  * Scargle (1982) Scargle J. D., 1982, ApJ, 263, 835
  * Schwope & Beuermann (1990) Schwope A. D., Beuermann K., 1990, A&A, 238, 173
  * Schwope et al. (1993) Schwope A. D., Beuermann K., Jordan S., Thomas H. C., 1993, A&A, 278, 487
  * Silber (1992) Silber A. D., 1992, PhD thesis, Massachusetts Institute of Technology
  * Sirotkin & Kim (2009) Sirotkin F. V., Kim W.-T., 2009, ApJ, 698, 715
  * Spruit & Ritter (1983) Spruit H. C., Ritter H., 1983, A&A, 124, 267
  * Thorstensen & Skinner (2012) Thorstensen J. R., Skinner J. N., 2012, The Astronomical Journal, 144, 81
  * Thorstensen et al. (2017) Thorstensen J. R., Ringwald F. A., Taylor C. J., Sheets H. A., Peters C. S., Skinner J. N., Alper E. H., Weil K. E., 2017, Research Notes of the AAS, 1, 29
  * Tovmassian et al. (1997) Tovmassian G. H., Greiner J., Zickgraf F. J., Kroll P., Krautter J., Thiering I., Zharykov S. V., Serrano A., 1997, A&A, 328, 571
  * Tremblay & Bergeron (2009) Tremblay P. E., Bergeron P., 2009, ApJ, 696, 1755
  * Verbunt (1984) Verbunt F., 1984, MNRAS, 209, 227
  * Verbunt & Zwaan (1981) Verbunt F., Zwaan C., 1981, A&A, 100, L7
  * Warner (1995) Warner B., 1995, Cataclysmic Variable Stars. Cambridge Astrophysics, Cambridge University Press, doi:10.1017/CBO9780511586491
  * Wright et al. (2010) Wright E. L., et al., 2010, AJ, 140, 1868
  * van Paradijs (1986) van Paradijs J., 1986, MNRAS, 218, 31p