# Revisiting the Conceptualization of Multiple Linear Regression

Grayson L. Baird
Brown University and Rhode Island Hospital
and
Stephen L. Bieber
University of Wyoming

###### Abstract

The problem known as multicolinearity has long been recognized to
fundamentally and negatively influence multiple regression. This paper does
not intend to either propose a numerical assessment of the degree to which
this problem exists within any data set or a solution to the problem itself.
Rather, it is our intent to illustrate the potentially serious ramifications
multicolinearity has on the traditional development of the multiple linear
regression (MLR) model and its associated statistics using established
equations, Venn diagrams, and real data.

Keywords: multicolinearity, collinearity, residualized regression, ordered
variable regression

##  1 Background

The problem known as multicolinearity has long been recognized to
fundamentally and negatively influence multiple regression. It is not the
intent of this paper to propose either a numerical assessment of the degree to
which this problem exists within any dataset or to propose a solution to the
problem itself. Rather, it is our intent to illustrate the potentially serious
ramifications multicolinearity has on the traditional development of the
multiple linear regression (MLR) model and its associated statistics.

The heart of the problem in the traditional development of MLR seems to stem
from the basic analysis of variance identity, in which the variation within
the dependent variable (Y) can be expressed as SS(total) = SS(model) +
SS(residual). Although this logically leads in regression to

$SS(total)=SS(regression)+SS(residual)$ (1)

, as seen in Woolf (1951), Baird and Bieber (2016, 2020), and below, equation
(1) does not hold for multiple regression problems in which the predictors are
correlated. The resultant impact of this violation can be rather extensive and
serious. As an aside, as in analysis of variance situations, equation (1) does
hold in regression contexts in which some orthogonalization of the variable
set is used (Yates, 1938). However, for any MLR model in which any or all of
the predictors (e.g., X1, X2) are correlated with each other and Y (Figure 1),
it will never hold. The calculation of the sum of squares (regression) =
SS(regression) is usually determined either directly as

$SS(regression)=\sum_{i=1}^{p}{\textup{b}_{i}}\sum_{i=1}^{p}{\textup{X}_{i}{Y}}={\textup{b}^{\prime}{X}^{\prime}{Y}}$
(2)

or indirectly once SS(residual) has been calculated as

$SS(regression)=\\\ SS(total)-SS(residual)$ (3)

. From this point onward, to avoid confusion, we will only refer to
SS(regression) and SS(residual) in association with their calculation through
equations (1), (2), and (3). The sum of squares associated with a particular
regression model will be symbolically referred to as SS(variables included).
For instance, the model sum of squares for the MLR model using the variables
X1 and X2 would be SS(X1,X2).

To facilitate the presentation of the basic principles and for the
confirmation of the assertions above, the Dwaine Studios, Inc. data used to
present multiple regression in Applied Linear Statistical Models (Kutner,
Nachtsheim, Neter, and Li; 2005, 5th edition, p. 237) will be used here; these
data and corresponding R code (R Foundation for Statistical Computing, Vienna,
Austria) are presented in Appendix A. However, any dataset and the development
of the fundamental statistics associated with MLR from any source would work
equally well for this discussion. The dependent variable (Sales) is
represented by Y throughout our discussion, and similarly the independent
variables (TARGTPOP and DISPOINC) will be represented by X1 and X2. Table 1
presents selected results from the simple regressions and the multiple
regression (MLR) using X1 and X2 to predict Y.

Figure 1: Relationship between Predictors (X1), (X2), and Dependent Variable
(Y) that are intercorrelated

Using the MLR coefficients presented in Table 1 along with the sums of squares
and cross-products provided in Table 2, the elements of equation (1) are

$SS(total)=\\\ Y^{\prime}Y=26,196.21$

$\displaystyle SS(regression)=b^{\prime}X^{\prime}Y=\
{\textup{b}_{1}}^{\prime}{\textup{X}_{1}}^{\prime}Y+{\textup{b}_{2}}^{\prime}{\textup{X}_{2}}^{\prime}Y$
$\displaystyle=(1.452)(12,730.59)+(9.366)(587.04)=24,015.28$

$\displaystyle
SS(residual)=SS(total)-SS(regression)=26,196.21-24,015.28=2,180.93$

These results above can be confirmed in Kutner, Nachtsheim, Neter, and Li
(2005) or by application of any statistical computer package. They provide the
foundation as a fait accompli for the presentation of either equation (2) or
(3) as the correct determination of the sum of squares associated with the MLR
model. The implication is that SS(regression) = SS(X1,X2) in the MLR
situation.

The question we would like to raise at this point is the following: is the
information provided above for equations (1), (2), and (3) actually correct?
More specifically, is the calculation of SS(regression) correct, in which it
is assumed that SS(regression) = SS(X1, X2)? In an attempt to answer this
question, we would like to digress for a moment and briefly introduce the work
presented in Woolf (1951). Woolf’s work is not being presented as the answer
or even as an answer to the multicolinearity problem, his work is merely being
used to assist us in answering the question above.

In 1951 computers were not readily available to researchers to make the
relatively complex calculations associated with multiple regression. To assist
researchers, Woolf (1951) wrote his paper on the Computation and
Interpretation of Multiple Regressions. His second method of calculating the
regression coefficients and analysis of variance sum of squares is known often
as the determination of Extra Sum of Squares.

Very simply, for our two predictor problem, the sum of squares for regression
is calculated as

$SS(regression)=SS(X\textsubscript{1})+SS(X\textsubscript{2}|X\textsubscript{1})$
(4)

where the predictors are ordered in their consideration. In equation (4), X1
is considered first and then X2. The Extra Sum of Squares associated with X2,
the sum of squares added to the explanation of Y beyond those already
established by X1, is calculated as the sum of squares associated with the
partial predictor (X2 |X1) as SS(X2|X1) =
SS(regression)–SS(X1).222SS(regression) as calculated by equation (2) or (3)
is intentional indicated here rather than SS(X1,X2). This concept is presented
in analysis of variance as the sequential determination of the sum of squares
known as Type I. From the information provided in Table 1, SS(X2 |X1) =
24,015.28 – 23,371.81 = 643.99. It is worth noting that:

$SS(regression)=SS(X\textsubscript{1})+SS(X\textsubscript{2}|X\textsubscript{1})=23,371.81+643.81=24,015.28$
(5)

This result is depicted in Figure 2 where regions 1 and 4 are accounted for by
SS(X1) and region 2 is accounted for by SS(X2 |X1).

Figure 2: Variation of Y Explained by the Orthogonal Predictors (X1) and (X2
|X1)

Not too surprisingly, these Type I sum of squares, when sequentially added,
will always equal SS(regression) as calculated by equation (2) or (3), since
these additions are always orthogonal. By changing the order in which he
placed the predictors in his second computational method, Woolf was able to
calculate the “fractions of variance associated with” (1951, p. 113) each of
the variables in the final MLR model. In effect, he showed that the sum of
squares associated with any predictor in the MLR model was SS(that predictor
|all the other predictors in the model), which is perhaps the first
recognizable presentation of the simultaneous estimation of the regression sum
of squares, Type III, for the regression application. Thus, for a two-
predictor situation such as the data in the Appendix, the sum of squares
associated with the MLR model is

$SS({\textup{X}_{1}},{\textup{X}_{2})}=SS({\textup{X}_{1}}|{\textup{X}_{2})}+SS({\textup{X}_{2}}|{\textup{X}_{1})}$
(6)

Now this result is surprising. The sum of the sequential orthogonal
predictors, equation (4), cannot possibly equal the sum of the simultaneous
partial predictors, equation (6), in any situation except when the predictors
are uncorrelated. This surprising result was also noted by Woolf, “Until now
analysis of variance (associated with the regression model)333The
parenthetical addition is included by the authors. has given sensible results.
But pushed too far, the method soon leads to apparent absurdity” (1951, p.
112). As a consequence, application of equation (1) for Woolf’s data resulted
in only 51 % of the SS(total) actually being accounted for. His conclusion was
decidedly disconcerting in 1951 and unfortunately is still disconcerting
today. As Woolf concluded,

> “There seems to be an enormous deficit in the balance sheet. The paradox, of
> course, is easily resolved. It is essential for the validity of an analysis
> of variance that the various components shall be uncorrelated or
> independent.” (1951, p. 113)

The actual sum of squares associated with the MLR model, SS(X1,X2), can be
calculated in several different ways. A few are presented below. The easiest
might be to simply calculate the Extra Sum of Squares contributions from
different starting points, such as Woolf did. For our two predictor situation,
using the data in Table 1 produces

$SS({\textup{X}_{2}}|{\textup{X}_{1})}=SS(regression)-SS({\textup{X}_{1}})=24,015.28-18,299.78=5,715.51$
(7)

$SS({\textup{X}_{1}}|{\textup{X}_{2})}=SS(regression)-SS({\textup{X}_{2}})=24,015.28-23,371.81=643.48$
(8)

and thus the actual sum of squares associated with the MLR model is

$SS({\textup{X}_{1}},{\textup{X}_{2})}=SS({\textup{X}_{1}}|{\textup{X}_{2})}+SS({\textup{X}_{2}}|{\textup{X}_{1})}=5,715.51+643.48=6,358.99$
(9)

Hence, the calculation of SS (regression) as shown on page 2 and in Table 1
does not equal SS(X1, X2), and the answer to the question we posed is no. The
sum of squares for the actual MLR model does not equal those produced by the
common calculation (i.e., equations 1, 2, 3). This result is perhaps better
explained through our second method of calculating the exact sum of squares
associated with the MLR model.

The second method is to first calculate the real data being used in the
multiple regression. The coefficients in the MLR model represent the
contribution of a particular predictor holding the other predictors constant.
In other words, the real data associated with the first predictor is not X1,
but rather X1|X2 = X1.2. These values are easily seen as the residuals from
the regression

${\textup{\^{X}}_{1}}=constant+{\textup{{a}}_{2}}{\textup{{X}}_{2}}$ (10)

and the values for X2 |X1 = X2.1 are the residuals from

${\textup{\^{X}}_{2}}=constant+{\textup{{a}}_{1}}{\textup{{X}}_{1}}$ (11)

. If these are the correct values, then the simple regressions of these
residualized predictors should reproduce the Extra Sum of Squares presented in
equations (7) and (8), and should reproduce the MLR regression coefficients
presented in Table 1. These simple regressions appear in Table 3 and confirm
these expectations. The result of this reassessment of the actual sum of
squares associated with the MLR model is that for the data of Table 1, the

$SS(\text{actual
total})=SS({\textup{X}_{1}},{\textup{X}_{2})}+SS(residual)=6,358.99+2,180.93=8,839.92$
(12)

This leaves 26,196.21 – 8,839.92 = 17,352.29 or 66.2 % missing from the total
sum of squares!

As seen in Figure 3, only regions associated with the residualized predictors,
regions 1 and 2, are included in Y. Region 4 has been entirely removed from
X1, X2, and Y. To quote Woolf, “There seems to be an enormous deficit in the
balance sheet.” (1951, p. 113). Unfortunately, unlike the pooling protocol in
analysis of variance, MLR does not place the interactive component, region 4,
in the residual when it is not included in the model. This confirms the
results in Figure 3: only the regions associated with the residualized
predictors and the residual are included in Y. Perhaps even more disconcerting
than the finding that the sum of squares associated with the MLR model is
incorrectly determined, is this finding that the application of the MLR model
does not reproduce the total variation of Y.

Figure 3: Variation of Y Explained by the Partial Predictors (X1.2) and (X2.1)

Are there other ramifications associated with the classical development of the
MLR model? For instance, even though the simple regressions of Table 3 have
replicated the unstandardized MLR coefficients, they have failed to produce
the same standardized coefficients. Which, if either, are correct? The
detailed calculation and explanation of the findings below can be found in
(Baird and Bieber, 2020) and hence only the conclusions are illustrated here.
To begin, there are two basic principles in regression that will hold for all
regression models in which equation (1) is valid. The first is

$\sum_{i=1}^{p}{\textup{z}_{i}^{2}}={\textup{R}^{2}}$ (13)

where zi is the ith standardized MLR coefficient. Several authors, such as
Darlington (1968), have recognized that equation (13) holds for data
situations in which the predictor variables are uncorrelated. Unfortunately, a
simplistic, additive, alternative to equation (13) has not been found for the
situation in which the predictors are correlated. However, rather than being
the result that the predictors are correlated, perhaps this inability has been
the result of using the incorrect data?

For the MLR model of the data in Table 1, if the sum of squares associated
with the model is wrong, then what is the correct R2? This is easily
calculated, within the context of the original problem, as SS(X1,X2)/SS(total)
= 6,358.99/26,196.21 = .243. Applying equation (13) to the standardized
coefficients in Table 1 produces (.748)2 \+ (.251)2 = .623, which is neither
the actual R2 nor the one presented in Table 1. However, applying equation
(13) to the standardized coefficients in Table 3 produces (.467)2 \+ (.157)2 =
.243, the correct R2. Therefore, the standardized coefficients in Table 3 are
logically correct. The use of the correct data which produces the correct
standardized coefficients validates the use of equation (13) for all MLR
situations regardless of whether a relationship exists within the predictor
variables or not. The second principle is the following:

$\sum_{i=1}^{p}{\textup{t}_{i}^{2}}/p={f}$ (14)

. Similarly, for the MLR model of the data in Table 1, what is the correct f?
Once again, this is easily calculated as MS(X1,X2) / MS(residual) =
[(6,358.99/2)/(2,180.93/18)] = 3,179.50/121.16 = 26.64, is the correct value,
not the value for f presented in Table 1. Applying equation (14) to the t
values presented in Table 1 produces [(6.87)2 \+ (2.31)2]/2 = 26.24, the
correct value (see Baird and Bieber, 2016; 2020). Thus, it can be seen that
the application of equation (2) in the MLR situation produces the incorrect
estimate of the regression sum of squares since it, as seen above, used the
wrong data. The residualized predictors, X1.2 and X2.1, from equations (10)
and (11), should be used rather than the original predictors X1 and X2. If
equation (2) is modified as the sum of the coefficients times the sum of the
cross-products between the residualized predictors and Y, then the correct
result is obtained. Using the unstandardized coefficients from Tables 2 and 4
and the Sum of Cross-Products from Table 2 gives

$\
{\textup{b}_{1}}\sum{\textup{X}_{1.2}}Y+{\textup{b}_{2}}\sum{\textup{X}_{2.1}}Y=(1.455)(3,929.37)+(9.366)(68.71)=6,358.99$
(15)

the value given in equation (9).

It can now be seen that the appropriate model associated with the Extra Sum of
Squares presentation use the orthogonal predictors X1 and X2.1, or X2 and
X1.2, rather than the original data. This can be seen in Table 4. These
orthogonal functions replicate perfectly the Extra Sum of Squares presentation
and are in fact the orthogonal functions presented in Woolf as his solution to
the paradox. They are his “formulae for the residuals.” (1951, p. 114) The
application of equation (2) to calculate SS(regression) is validated by either
of the two orthogonal regressions in Table 4.

In conclusion, Table 1 provides the traditional summary for an MLR model. As
seen above, when the predictors are correlated, only the unstandardized
coefficients, the t-values, and the SS(residual) of this table are correct.
Unfortunately, the SS(regression), the R2, the f, and the standardized
coefficients are all incorrect. It is also unfortunate that the discussion of
the Extra Sum of Squares, Type I determination of component sum of squares, is
linked inappropriately with the Type III MLR model. As a consequence of these
errors, both SPSS (IBM SPSS Statistics, Armonk, NY) and SAS (SAS Software,
Cary, NC) use the traditional multiple regression presentation to develop
their analysis of variance tables, identical to our Table 1. In an interesting
twist, base R (R Foundation for Statistical Computing, Vienna, Austria)
follows the Extra Sum of Squares presentation and provides a Type I analysis
of variance table, but the parameter estimates from the MLR Type III model.

##  2 References

  1. 1.

Baird, G. L., & Bieber, S. L. (2016). The Goldilocks dilemma: Impacts of
multicollinearity-a comparison of simple linear regression, multiple
regression, and ordered variable regression models. Journal of Modern Applied
Statistical Methods, 15(1), 18.

  2. 2.

Baird, G. L., & Bieber, S. L. (2020). Sampling the Porridge: A Comparison of
Ordered Variable Regression with F and R 2 and Multiple Linear Regression with
Corrected F and R 2 in the Presence of Multicollinearity. Journal of Modern
Applied Statistical Methods, 18(1), 11.

  3. 3.

Darlington, R. B. (1968). Multiple regression in psychological research and
practice. Psychological bulletin, 69(3), 161.

  4. 4.

Fox, J., Weisberg, S., Adler, D., Bates, D., Baud-Bovy, G., Ellison, S.,
Firth, D., Friendly, M., Gorjanc, G., Graves, S. & Heiberger, R., Package
‘car’. (2022). Version R package version 4.1.2.
https://cran.r-project.org/web/packages/car/index.html

  5. 5.

French, J., Package ‘api2lm’. (2022). Version R package version 4.1.2.
https://cran.r-project.org/web/packages/api2lm/api2lm.pdf

  6. 6.

Kutner, M. H., Nachtsheim, C. J., Neter, J., & Li, W. (2005).Applied Linear
Statistical Models, McGraw Hill Irwin, New York. NY, 409.

  7. 7.

Revelle, W., & Revelle, M. W. (2015). Package ‘psych’. The comprehensive R
archive network, 337, 338.

  8. 8.

Woolf, B. (1951). Computation and interpretation of multiple regressions.
Journal of the Royal Statistical Society: Series B (Methodological), 13(1),
100-119.

  9. 9.

Yates, F. (1938). Orthogonal functions and tests of significance in the
analysis of variance. Supplement to the Journal of the Royal Statistical
Society, 5(2), 177-180.

Table 1: Simple and Multiple Linear Regression Results for Y, X1 , X2

| Simple (X1) |  |  |  | Simple (X2) |  |  |  | MLR (X1,X2) |  |
---|---|---|---|---|---|---|---|---|---|---|---
Source | SS | F | R2 |  | SS | F | R2 |  | SS | F | R2
Reg | 23,371.81 | 157.22 | .892 |  | 18,299.78 | 44.03 | .699 |  | 24,015.28 | 99.10 | .917
Res | 2,824.40 |  |  |  | 7,896.43 |  |  |  | 2,180.93 |  |
Total | 26,196.21 |  |  |  | 26,196.21 |  |  |  | 26,196.21 |  |
| b | z | t |  | b | z | t |  | b | z | t
X1 | 1.836 | .945 | 12.54 |  |  |  |  |  | 1.455 | .748 | 6.87
X2 |  |  |  |  | 31.173 | .836 | 6.64 |  | 9.366 | .251 | 2.31

  * •

Note: z is the standardized coefficient

Table 2: Sum of Squares and Cross-Products
(Mean Centered) | Y | X1 | X2
---|---|---|---
Y | 26,196.21 |  |
X1 | 12,730.59 | 6,934.33 |
X2 | 587.04 | 282.33 | 18.83
X1.2 | 3,929.37 |  |
X2.1 | 68.71 |  |
Table 3: Simple Regression Results for Y, X1.2, X2.1

| Simple (X1.2) |  |  |  | Simple (X2.1) |  |
---|---|---|---|---|---|---|---
Source | SS | f | R2 |  | SS | f | R2
Regression | 5,715.51 | 5.30 | .218 |  | 643.48 | .48 | .025
Residual | 20,480.71 |  |  |  | 25,552.73 |  |
Total | 26,196.21 |  |  |  | 26,196.21 |  |
| b | z | t |  | b | z | t
X1.2 | 1.455 | .467 | 2.30 |  |  |  |
X2.1 |  |  |  |  | 9.366 | .157 | .69

  * •

Note: z is the standardized coefficient

Table 4: Orthogonal Function Regression Results: Y, X1, X2.1 and Y, X2, X1.2

| X1 and X2.1 |  |  |  | X2 and X1.2 |  |
---|---|---|---|---|---|---|---
Source | SS | f | R2 |  | SS | f | R2
Regression | 24,015.28 | 99.10 | .917 |  | 24,015.28 | 99.10 | .917
Residual | 2,180.93 |  |  |  | 2,180.93 |  |
Total | 26,196.21 |  |  |  | 26,196.21 |  |
| b | z | t |  | b | z | t
X1 | 1.836 | .945 | 13.89 |  |  |  |
X2.1 | 9.366 | .157 | 2.31 |  |  |  |
X2 |  |  |  |  | 31.173 | .836 | 12.29
X1.2 |  |  |  |  | 1.455 | .467 | 6.87

  * •

Note: z is the standardized coefficient

##  3 Appendix A

    #R Code

    #Dataset from Kutner, et al. (2005) :
    library("api2lm")

    X1=dwaine$targetpop
    X2=dwaine$dpi
    Y= dwaine$sales
    ###Table 1

    ## SLR Models
    SLR1 = lm(Y~ X1)

    #ANOVA
    anova(SLR1)

    #Coefficients
    summary(SLR1)
    B1= (1.8359*(sd(X1))/sd(Y))
    print(B1)

    SLR2 = lm(Y~ X2)

    #ANOVA
    anova(SLR2)

    #Coefficients
    summary(SLR2)
    B2= (31.173*(sd(X2))/sd(Y))
    print(B2)

    ### MLR Model
    MLR = lm(Y~ X1 + X2)

    #ANOVA
    library("car")
    Anova(MLR, type="III")
    anova(MLR) #type I

    #Coefficients
    summary(MLR)
    B1= (1.4546*(sd(X1))/sd(Y))
    print(B1)
    B2= (9.3655*(sd(X2))/sd(Y))
    print(B2)

    ###Table 3

    ##Residualized Predictors
    residX1 = lm(X1~X2)
    X1.2=resid(residX1)

    #SLR Models
    SLR1.2 = lm(Y~ X1.2)
    summary(SLR1.2)
    anova(SLR1.2)

    residX2 = lm(X2~X1)
    X2.1=resid(residX2)

    SLR2.1 = lm(Y~ X2.1)
    summary(SLR2.1)
    anova(SLR2.1)

    #Correct Standardized Coefficients and Correct R2
    B1.2= (1.4546*(sd(X1.2))/sd(Y))
    print(B1.2)
    Correct_R2=(B1.2^2)
    print(Correct_R2)

    B2.1= (9.3655*(sd(X2.1))/sd(Y))
    print(B2.1)
    Correct_R2=(B2.1^2)
    print(Correct_R2)

    ###Table 2 SSCP
    library("psych")
    all=data.frame(Y,X1,X2,X1.2,X2.1)
    sscp=cov(all)/.05
    print(sscp)

    ###Table 4
    ##Orthogonal Functions Model
    OFM1 = lm(Y~ X1+X2.1)
    anova(OFM1)

    #Coefficients and R2
    summary(OFM1)

    B1= (1.8359*(sd(X1))/sd(Y))
    print(B1)

    B2.1= (9.3655*(sd(X2.1))/sd(Y))
    print(B2.1)

    R2=(B1^2)+(B2.1^2)
    print(R2)

    ###Orthogonal Functions Model 2
    OFM2 = lm(Y~ X2+X1.2)
    anova(OFM2)

    ##Coefficients ND R2
    summary(OFM2)

    B2= (31.1732*(sd(X2))/sd(Y))
    print(B2)

    B1.2= (1.4546*(sd(X1.2))/sd(Y))
    print(B1.2)

    R2=(B2^2)+(B1.2^2)
    print(R2)

    ###Equation 12
    MLR = lm(Y~ X1+X2)

    #Type III Sums of Squares
    Anova(MLR, type="III")
    Actual_Total_SS=(5715.5+643.5+2180.9)
    print(Actual_Total_SS)

    ###Equation 13
    B1..2= (1.4546*(sd(X1.2))/sd(Y))
    B2..1= (9.3655*(sd(X2.1))/sd(Y))
    Correct_R2=(B1..2^2)+(B2..1^2)
    print(Correct_R2)


# Detection Threshold of Audio Haptic Asynchrony in a Driving Context

Gyanendra Sharma
Toyota Research Institute
Los Altos, CA 94022
<EMAIL_ADDRESS>
Hiroshi Yasuda
Toyota Research Institute
Los Altos, CA 94022
<EMAIL_ADDRESS>
Manuel Kuehner
Toyota Research Institute
Los Altos, CA 94022
<EMAIL_ADDRESS>
###### Abstract

In order to provide perceptually accurate multimodal feedback during driving
situations, it is vital to understand the threshold at which drivers are able
to recognize asyncrony between multiple incoming Stimuli. In this work, we
investigated and report the detection threshold (DT) of asynchrony between
audio and haptic feedback, in the context of a force feedback steering wheel.
We designed the experiment to loosely resemble a driving situation where the
haptic feedback was provided through a steering wheel (Sensodrive), while the
accompanying audio was played through noise cancelling headphones. Both
feedbacks were designed to resemble rumble strips, that are generally
installed on the side of major roadways as a safety tool. The results indicate
that, for $50\%$ of the participants, asynchrony was detectable outside the
range of -75 ms and 110 ms, where the former is related to perceiving audio
before haptic and vice versa for the latter. We were also able to concur with
previous studies, which state that latency is perceivable at a lower threshold
when audio precedes haptic stimuli.

##  1 Introduction

Rumble strips on the side of the road are a major safety tool to make sure
drivers do not drift off the road during driving situations. These physical
artifacts assert driver’s attention by providing both auditory and haptic
feedback. Designing safety systems that translate similar phenomena in a
synthetic manner, so that drivers can be warned in cases of road departures
holds significant promise. Specifically, an interaction concept that offers
digitally created rumble strip phenomena through audio and haptic feedback can
be created to warn drivers when they unintentionally depart the lane or the
road itself. One key advantage of such a concept is an easy translation of
mental model from an already existing warning mechanism that is widely
understood and familiar, especially in North America.

Figure 1: Rumble Strip in a physical environment.

In order to do this, there are various parameters that characterize the
feedback, and one of the major components is the timing at which both audio
and haptic responses are delivered to the driver. While physical rumble
strips, as shown in figure 1 are characterized by natural law of physics,
synthetic or electronically created feedback systems require each signal to
propagate through a separate set of actuators, control units and computation
units before they arrive to the driver. This causes delay in regards to the
delivery of either of the haptic or audio stimuli to the driver. Figure 2
shows the case where haptic stimuli is delivered to the driver slightly later
than audio. A large asynchrony between the stimuli representing the same
feedback event can lead to confusing user feedback systems. So, for these
feedback systems to be effective and also to avoid annoyance, which has been
shown to be a major reason for safety systems being turned off, it is
pertinent that both audio and haptic signal onset time are perceptually in
sync, even if asynchrony exists in the underlying system. The detection
threshold (DT) at which this asynchrony is perceivable to the driver is the
major contribution of this work.

Figure 2: Latency between audio and haptic stimuli in the context of real
world vs synthetic rumble strips.

Human beings can only detect asynchrony outside of a certain time window. This
is also referred to as the point of subjective simultaneity (PSS), where
stimulus onset asynchrony (SOA) are most likely to yield responses that affirm
that the SOA was detectable Stone et al. (2001). Early works on this topic are
derived from the field of psychology, where SOA between audio and visual
stimuli were studied Hershenson (1962); Woodworth et al. (1954). More
recently, haptic has also been included to study the perception of cross-modal
simultaneity Adelstein et al. (2003); Levitin et al. (2000); Altinsoy (2003).
Beyond the reported DT numbers, these studies also demonstrated that latency
is perceived at a lower threshold when audio precedes the haptic feedback.

While earlier studies focused on a more physical setup, recent studies have
investigates this topic in the context of modern applications such as specific
computer interfaces, touchscreens and force feedback haptic devices. For
instance, work of Silva et al. Silva et al. (2013) was specifically focused on
computer based multimedia application and showed that asynchrony of audio
played before haptic within 92 ms and, audio played after haptic within 110 ms
was not detected by the participants. With the rapid increase in usage of AR
and VR technologies in recent years, the application areas upon which this
topic is being investigated has expanded significantly. Tele-operation
Neumeier et al. (2019); Liu et al. (2017), multi-sensory media (mulsemedia)
Covaci et al. (2018); Yuan et al. (2015) are a few examples.

Through this user study, we investigate the detection threshold (DT) at which
drivers are able to perceive asynchrony between haptic and audio. The DT is
determined using a 3 alternative forced choice (3AFC) adaptive staircase
procedure Prins et al. (2016); Gescheider (1997). Threshold obtained in this
manner corresponds to the 50 percentile point on a psychometric function
Wichmann and Hill (2001).

In particular, we are interested in the following research questions.

  * •

What is the time difference at which drivers begin to notice asynchrony
between haptic and audio feedback.

  * •

Is the previous conclusion, latency is perceived at a lower threshold when
audio precedes the haptic, applicable to driving context?

##  2 System Characterization

Before embarking on designing a system that is able to accurately provide
desired offsets between audio and haptic feedback, we had to explore the role
of system latency i.e the inherent latency within the computing environment of
the workstation. We investigated this by recording repeated samples of audio-
haptic feedback pairs, as shown in figure 3, and noting the time difference
through analysis of their respective frequencies.

Figure 3: Setup to record audio and vibration.

We used a total of 58 audio-haptic pairs and analyzed each pair to determine
their actual offset time and compared it with the desired offset time to
obtain a value for system latency. The audio had variety of offsets ranging
between 6 ms to 30 ms and it was played at 1KHz while the vibration was played
at 100 Hz. As shown in figure 4, a spectrogram for each feedback pair was
analyzed to determine the actual offset. The audio onset time ($ao$), haptic
onset time ($ho$), desired offset ($do$) from all samples ($N=58$) was used to
calculate the system latency ($sl$) as follows.

Figure 4: Spectrogram to analyze the offset between audio and haptic.

$sl=\frac{1}{N}{\sum_{\alpha=1}^{N}((ao_{\alpha}-ho_{\alpha})-do_{\alpha})}$
(1)

Using equation 1, we calculated the system latency on average to be +15.93 ms,
with median value of 16 ms and standard deviation of 4.61 ms. This showed that
audio on average gets delayed by approximately 15.93 ms due to system latency
even in cases when both are played simultaneously. We included the finding in
our experimental design and setup to compensate for this system delay.

##  3 Methodology

In this study, participants held onto a Sensodrive steering wheel sen (2022)
that was connected to a Linux workstation while wearing a pair of noise
canceling wired headphones. Participants were exposed to 3 sets of haptic and
audio feedback pairs with 3 seconds interval between each set. The haptic
feedback appeared in the form of vibration on the Sensodrive wheel while the
audio feedback was played through the headphone. Two of the audio-haptic
feedback pairs were in sync while the third one was played out of sync with
experimentally controlled offset values. The order at which each of these sets
appeared was randomized. The pattern of the feedback for both the sound and
haptic were designed to resemble a rumble strip. After each trial,
participants were asked the following question to understand their perception
of the feedback; “which audio-haptic pair do you think was NOT in sync?”

We employed a 3AFC (3 alternatives force choice) method i.e participants can
respond with only one of the three given choices. In this particular case,
they were instructed to respond by speaking out one of the three numbers; 1, 2
or 3 based on whether the first, second or third feedback pair was out of
sync. Their responses were recorded by the experimenter and the next trial
would commence. The exposure to stimulus followed the 2 step up, 1 step down
staircase approach as shown in figure 5. So, the stimulus separation between
audio and haptic for the out of sync feedback pair is determined based on
user’s previous input. For successive correct responses, the stimulus
difference would be decreased, whereas a single incorrect response would lead
to increase in stimulus latency between the audio haptic pair. This leads to
each participant responding in ways where they move back and forth within a
certain time window of stimulus separation. Each participant was continuously
exposed in this manner until 8 reversal of responses was recorded. A
representation of an actual participant response is shown in 5.

Figure 5: Sample participant response based on the 2 step up, 1 step down
staircase approach.

We designed a within subjects study since we wanted to expose each participant
to both latency conditions; audio preceding haptic and vice versa. Each
participant underwent two experimental conditions; (i) negative audio haptic
latency, where audio is played before the haptic feedback and (ii) positive
audio haptic latency, where audio is played after the haptic feedback. To
overcome ordering influence, each participant was randomly assigned to one of
the experimental conditions to begin with and moved to the next one. Once the
participants completed this stage, they were asked to complete a post task
questionnaire to gather information on driving experiences, gender, age and
other possible relevant information.

###  3.1 Pilot Study

We conducted pilot studies with participants that were employees of an
automotive company. We used 4 participants to explore and establish certain
parameters that were later applied to the user study. This included the
starting stimulus offset for both haptic or audio first setup, number of
reversal points, and the step size. We initially overestimated the starting
stimulus offset for both audio and haptic first setup to be between 350 ms -
400 ms. However, it became quite obvious during the pilot studies that a more
conservative estimation of approximately 260 ms was adequate. In addition,
during the pilot studies, we used 9 reversal points to end the experiment for
each session. However, we realized that, based on the length of the experiment
and overall convergence of stimulus values, decreasing it to 8 reversals,
excluding the first reversal in both cases, was more appropriate. In terms of
step size, our initial estimate of 30ms as large or initial step size and 15
ms as the small step size didn’t seem to cause any issues and therefore, we
continued with the same numbers for the user study.

###  3.2 Participants

We recruited 15 external participants to conduct the user study out of which 7
were female and 8 were male. In terms of the age distribution, 7 were above
the age of 40 and 8 below. In addition, all participants were experienced
drivers, with 11 participants indicating that they had more than 10 years of
driving experience, with only 2 participants indicating less than 5 years of
driving experience. 3 participants indicated that they drove 4-6 times a week
while the rest indicated that they drove daily.

###  3.3 Experimental Procedure

Figure 6: Setup with the force feedback sensodrive wheel and headphone
connected to the workstation. Not seen in the figure - Linux workstation.

All participants were hired using the services of Fieldwork San Francisco. The
entire experiment was conducted over a period of approximately one week. Once
each participant arrived at the research facility, they went through all the
covid protocols put in place before they were allowed inside. Each participant
signed the consent form and also verbally consented to being recorded before
they were given a brief description about the experiment and their task.
Before commencing the experiment, each participant went through a simple test
to make sure the sound being played on the wired headset was at an acceptable
level while also making sure that they were able to feel the vibration on the
steering wheel. The physical setup is shown in figure 6.

As a safety protocol, we instructed the participants to hold the steering
wheel in a 10 - 2 position and make sure to not have their thumb or fingers
inside the wheel to avoid injuries in an off chance that the wheel makes
sudden movements. At this time, participants were allowed to ask any questions
they might have before beginning the trials. After a session was concluded i.e
8 recorded reversals, participant were told that they could rest their hands
for a minute before beginning the second session. Once this was over, they
were instructed to fill out a post task questionnaire and this marked the end
of the experimental session. Any thoughts, questions or concerns they had were
verbally discussed and noted by the experimenter at this time before escorting
them out of the research facility.

##  4 Results

Figure 7: Results of stimulus responses for all 15 participants in cases where
haptic precedes audio. Figure 8: Results of stimulus responses for all 15
participants in cases where audio precedes haptic. Figure 9: Percentage of the
participants who could detect the asynchrony. The red dotted line indicates
the detection threshold for $50\%$ of participants.

The results for each participant is shown in figure 7 and 8. In terms of the
analysis of the overall results, participants 1 and 3 were discarded as
outliers as their responses did not follow even the minimal expected pattern
i.e they consistently failed to notice stimulus offsets that were larger than
even the first stimulus value in both the audio first as well as the haptic
first scenario. The DT was derived using equation 2 for the remaining 13
participants. When haptic feedback preceded audio, the average DT was reported
at $123.41\pm 4.61$ ms and $87.73\pm 4.61$ ms when audio preceded haptic. We
were unable to observe any effects of age and gender in terms of the results.

$DT=\frac{1}{7}{\sum_{\alpha=2}^{8}R_{\alpha}}$ (2)

The cumulative score of participants based on their individual DTs is shown in
figure 9. This shows that in $50\%$ of cases, the detection threshold lies
between -75 ms and 100 ms. Similarly, extrapolating from the same figure, we
can also see that $20\%$ of the participants/drivers would notice a latency of
50 ms when audio is played before haptic and 65 ms when audio is played after
haptic feedback. The implications based on these findings are briefly
discussed below.

##  5 Discussion

The results indicate that the detection threshold for a haptic and audio based
multi-modal feedback system in a driving context closely follow the findings
from the past i.e audio preceding haptic is observable at smaller absolute
time differences than the other way round.

On the other hand, there are also design implications that can be derived
based on the results reported above. For instance, safety systems in cars that
rely on providing audio-haptic feedback, the robustness of the system can be
tuned to accommodate the latency parameters based on our findings. In
situations where, if for example, we want to design a feedback system that
would only allow for $20\%$ of the participants to notice the asynchrony,
latency window between -50 ms to 65 ms could be used, based on extrapolation
on the graph shown in figure 9. Similarly, an HMI interface that requires no
more than 20% of the drivers to perceive the asynchrony will have to optimize
their computation power to include delivery of haptic and audio stimuli within
the range of -50 ms to 65 ms. Based on the computational cost as well as the
objective of the design, this can be tuned to either cover for more or less of
the population. Our results provide a clear guideline in terms of what these
human factors costs will amount to be.

##  6 Limitations and Future Work

This study was conducted entirely within a research setting and without the
use of real driving situations or scenarios. This is a major limitation in
regards to drawing accurate conclusions as the primary task for the
participants in these scenarios was not driving. So, the DT estimates derived
from this user study should be considered as conservative. We expect the DT
for the same set of feedback; haptic and audio to be significantly higher if a
similar setup were to be tested in a real driving context as cognitive
resources are expended on higher priority tasks. However, to what extent such
effects are characterized would be part of future research. The results
presented in this study can be applied as baselines to draw meaningful
extrapolations as needed.

The second major limitation was characterized by the Sensodrive wheel that was
used as the primary component to provide haptic feedback. While the audio
component is generated entirely in a digital format, that’s not exactly the
case with the Sensodrive wheel which used the physical motors to generate the
vibration. This in turn also creates sound artifacts that are clearly audible.
The usage of over the ear noise canceling headphones mitigate this factor but
does not entirely nullify the sound factor emanating from the wheel.

Overall, we were able to obtain results that we believe will have significant
implications in designing multi-modal feedback systems as part of vehicle
safety system framework. Our results, especially the one presented in figure
9, can be applied to characterize system specifications when steering based
multi-modal feedback systems are designed. Also, as evidenced by our results,
more care can be put in place to make sure audio isn’t played much earlier
than other modalities as threshold of latency perception in such scenarios is
lower.

## References

  * Stone et al. (2001) JV Stone, NM Hunkin, J Porrill, R Wood, V Keeler, M Beanland, M Port, and NR Porter.  When is now? perception of simultaneity.  _Proceedings of the Royal Society of London. Series B: Biological Sciences_ , 268(1462):31–38, 2001.
  * Hershenson (1962) Maurice Hershenson.  Reaction time as a measure of intersensory facilitation.  _Journal of experimental psychology_ , 63(3):289, 1962.
  * Woodworth et al. (1954) Robert Sessions Woodworth, Bernard Barber, and Harold Schlosberg.  _Experimental psychology_.  Oxford and IBH Publishing, 1954.
  * Adelstein et al. (2003) Bernard D Adelstein, Durand R Begault, Mark R Anderson, and Elizabeth M Wenzel.  Sensitivity to haptic-audio asynchrony.  In _Proceedings of the 5th international conference on Multimodal interfaces_ , pages 73–76, 2003.
  * Levitin et al. (2000) Daniel J Levitin, Karon MacLean, Max Mathews, Lonny Chu, and Eric Jensen.  The perception of cross-modal simultaneity (or “the greenwich observatory problem” revisited).  In _AIP Conference Proceedings_ , volume 517, pages 323–329. American Institute of Physics, 2000.
  * Altinsoy (2003) M Ercan Altinsoy.  Perceptual aspects of auditory-tactile asynchrony.  In _Proceedings of the Tenth International Congress on Sound and Vibration_ , pages 3831–3838. Citeseer, 2003.
  * Silva et al. (2013) Juan M Silva, Mauricio Orozco, Jongeun Cha, Abdulmotaleb El Saddik, and Emil M Petriu.  Human perception of haptic-to-video and haptic-to-audio skew in multimedia applications.  _ACM Transactions on Multimedia Computing, Communications, and Applications (TOMM)_ , 9(2):1–16, 2013.
  * Neumeier et al. (2019) Stefan Neumeier, Philipp Wintersberger, Anna-Katharina Frison, Armin Becher, Christian Facchi, and Andreas Riener.  Teleoperation: The holy grail to solve problems of automated driving? sure, but latency matters.  In _Proceedings of the 11th International Conference on Automotive User Interfaces and Interactive Vehicular Applications_ , pages 186–197, 2019.
  * Liu et al. (2017) Ruilin Liu, Daehan Kwak, Srinivas Devarakonda, Kostas Bekris, and Liviu Iftode.  Investigating remote driving over the lte network.  In _Proceedings of the 9th international conference on automotive user interfaces and interactive vehicular applications_ , pages 264–269, 2017.
  * Covaci et al. (2018) Alexandra Covaci, Longhao Zou, Irina Tal, Gabriel-Miro Muntean, and Gheorghita Ghinea.  Is multimedia multisensorial?-a review of mulsemedia systems.  _ACM Computing Surveys (CSUR)_ , 51(5):1–35, 2018\.
  * Yuan et al. (2015) Zhenhui Yuan, Ting Bi, Gabriel-Miro Muntean, and Gheorghita Ghinea.  Perceived synchronization of mulsemedia services.  _IEEE Transactions on Multimedia_ , 17(7):957–966, 2015.
  * Prins et al. (2016) Nicolaas Prins et al.  _Psychophysics: a practical introduction_.  Academic Press, 2016.
  * Gescheider (1997) GA Gescheider.  Psychophysics: The fundamentals. lawrence erlbaum associates.  _Inc., Publishers_ , pages 1–71, 1997.
  * Wichmann and Hill (2001) Felix A Wichmann and N Jeremy Hill.  The psychometric function: I. fitting, sampling, and goodness of fit.  _Perception & psychophysics_, 63(8):1293–1313, 2001.
  * sen (2022) Force feedbackc sensodrive wheel.  https://www.sensodrive.de/products/senso-wheel-sd-lc_e.php, 2022\.  [Online; accessed 02-Feb-2022].


# Circular Polarimetry of Suspect Wind-accreting Magnetic pre-Polars

Pasi Hakala1, Steven G. Parsons2, Thomas R. Marsh3, Boris T. Gänsicke3, Gavin
Ramsay4, Axel Schwope5, J.J. Hermes6
1Finnish Centre for Astronomy with ESO (FINCA), Quantum, University of Turku,
FI-20014, Finland
2Department of Physics and Astronomy, University of Sheffield, Sheffield, S3
7RH, UK
3Department of Physics, University of Warwick, Coventry CV4 7AL, UK
4Armagh Observatory and Planetarium, College Hill, Armagh, BT61 9DG, UK
5 Leibniz Institute for Astrophysics Potsdam (AIP), An der Sternwarte 16,
14482 Potsdam, Germany
6Department of Astronomy, Boston University, 725 Commonwealth Ave., Boston, MA
02215, USA
E-mail<EMAIL_ADDRESS>
###### Abstract

We present results from a circular polarimetric survey of candidate detached
magnetic white dwarf – M dwarf binaries obtained using the Nordic Optical
Telescope, La Palma. We obtained phase resolved spectropolarimetry and imaging
polarimetry of seven systems, five of which show clearly variable circular
polarisation. The data indicate that these targets have white dwarfs with
magnetic field strengths $>80$ MG. Our study reveals that cyclotron emission
can dominate the optical luminosity at wavelengths corresponding to the
cyclotron emission harmonics, even in systems where the white dwarfs are only
wind-accreting. This implies that a very significant fraction of the the
stellar wind of the companion star is captured by the magnetic white dwarf
reducing the magnetic braking in pre-CVs. Furthermore, the polarimetric
confirmation of several detached, wind-accreting magnetic systems provides
observational constraints on the models of magnetic CV evolution and white
dwarf magnetic field generation. We also find that the white dwarf magnetic
field configuration in at least two of these systems appears to be very
complex.

###### keywords:

Physical data and processes: accretion, magnetic accretion – stars: binaries
close – stars: magnetic fields

††pagerange: Circular Polarimetry of Suspect Wind-accreting Magnetic pre-
Polars–References

##  1 Introduction

Cataclysmic variables (CVs) are semi-detached binaries where a cool donor,
usually a main-sequence star (i.e. the secondary) fills its Roche-lobe and
matter is transferred via Lagrangian L1 point overflow and accreted by the
primary white dwarf (WD). In the case of non-magnetic WDs, an accretion disc
is formed, but if the WD magnetic field is strong enough ($\sim 1$ MG or
higher) the disc is at least partly disrupted and the inner accretion flow
onto the surface of the magnetic WD (MWD) follows the magnetic field lines. In
case of strong magnetic fields (i.e. $\sim 10$ MG or more), the entire
accretion disc is disrupted and the ballistic accretion stream leaving the L1
point will be eventually re-directed along the magnetic field lines of the MWD
(within its Alfén radius) and impacts the MWD near the magnetic pole(s) (see
Cropper 1990 for an early review). Another consequence of the strong magnetic
field of the WD is the synchronicity of the MWD spin period with the orbital
period of the system: these synchronous systems are called polars (or AM Her
binaries).

In addition to affecting the accretion geometry of these systems, the magnetic
fields of both the WD and the secondary star play a vital role in the binary
evolution. These sources initially form as a consequence of common envelope
(CE) phase of binary evolution (Paczyński, 1967) as one of the binary
components is swallowed by the outer envelope of the other, evolved, binary
component. This leads to a rapid loss of angular momentum and thus the binary
orbit spirals in. Finally the red giant phase ends, and as a result, a binary
consisting of a main sequence star and a compact object emerges. These systems
have periods from hours up to several days. The orbit of this binary shrinks
further as angular momentum is lost due to both magnetic braking and
gravitational radiation (see Paczyński 1967 and Verbunt & Zwaan 1981). At some
point the low mass main-sequence star will start to fill its binary Roche lobe
and, as the orbital shrinkage continues, mass will start to overflow the
Lagrangian L1 point and accrete on the primary component giving rise to a CV
in the case of a WD or an X-ray Binary for a neutron star or black hole.

However, in Low Accretion Rate Polars (LARP, Schwope et al. 2002) the
secondary star does not quite fill its Roche lobe, but there is still some
indication of matter being accreted on the WD as beamed cyclotron emission is
observed on/above the WD surface in the optical light curves. The hot spots
are indicative of a strong magnetic field of the WD, that channels the
accretion onto its magnetic poles. As there is no Roche lobe overflow, the
implication is that the accretion must take place via stellar wind from the
secondary star that is captured by the WD magnetosphere. These systems have
also later been called PREPs (i.e. pre-polars, Schwope et al. 2009) following
the realisation that most of them do not contain a Roche lobe filling donor.

Source | Date | Filter/grism | Mode | UT range | Exp. time (sec) | Notes
---|---|---|---|---|---|---
2MASS J01294349+6715301 | 16.10.2020 | Grism 5 | circ.sp.pol. | 20:46-06:31 | 180 | (three 20-30min gaps)
SDSS J222918.95+185340.22 | 17.10.2020 | Grism 5 | circ.sp.pol. | 19:37-19:46 | 200 | (4 circ. pol. spectra)
| 17.10.2020 | $R_{\mathrm{Bess}}$ | circ.im.pol. | 19:52- 23:14 | 30 |
| 16.11.2020 | $R_{\mathrm{Bess}}$ | circ.im.pol. | 20:43-22:35 | 20 |
SDSS J030856.55$-$005450.73 | 18.10.2020 | Grism 5 | circ.sp.pol. | 00:24-00:45 | 300 | (4 circ. pol. spectra, very faint)
| 18.10.2020 | $R_{\mathrm{Bess}}$ | circ.im.pol. | 00:50-05:05 | 30 | (some 60sec due to bad weather)
ZTF J014635.73+491443.14 | 11.12.2020 | Grism 10 | circ.sp.pol. | 23:21-01:48 | 500 |
| 15.12.2020 | Grism 10 | circ.sp.pol. | 23:48-02:25 | 180 |
SDSS J075015.11+494333.22 | 04.03.2021 | Grism 5 | circ.sp.pol. | 20:59-01:24 | 300 |
SDSS J085336.03+072033.52 | 04.03.2021 | $R_{\mathrm{Bess}}$ | circ.im.pol. | 20:03-00:05 | 30 |
SDSS J121209.31+013627.75 | 08.03.2021 | $i$ | circ.im.pol. | 01:42-02:34 | 45 | (cirrus)
Table 1: The Observing log where we indicate the source name; the start date
of the observations; the filter in the case of imaging polarimetry or Grism in
the case of spectroscopic polarimetry; the UT time range of the observations;
the exposure time of each image and any notes on the observations.(1Krushinsky
et al. (2020),2Parsons et al. (2021),3Becker et al. (2011),4Guidry et al.
(2021),5Schmidt et al. (2005b)

Recently Schreiber et al. (2021) performed detailed computations of magnetic
CV evolution, where they suggested that the stronger magnetic field strengths
observed in polars, as compared to single WDs, could be due to WD magnetic
field generation by a dynamo process (Isern et al., 2017), as the crystallized
WD is spun up by accreted matter. They managed to explain the observed
properties of magnetic CV population (and the lower field strengths observed
in single WDs) remarkably well. In their evolutionary model the detached pre-
polars are in fact a product of a _non-magnetic_ CV phase, preceding the polar
phase. These systems are predicted to populate certain orbital period
intervals and thus building up the observed population is crucial for testing
the model. This has also possible implications for the magnetic braking theory
of the binary period evolution (Li, Wu, & Wickramasinghe, 1994). If a large
fraction of the stellar wind is not carried out of the system, the magnetic
braking is at least greatly reduced, slowing down the evolution of the binary
towards shorter orbital periods and potentially making the observed population
of CVs much older.

For this study, we have compiled a sample of candidate "pre-polars", where the
secondary star does not fill its Roche lobe. The selected sources all have
previous observed properties suggesting the presence of possible cyclotron
emission (i.e. evidence for cyclotron humps in the optical spectrum;
extraordinary colours in light curves or excess emission in some passband).
The individual sources are introduced and discussed separately in the
subsections below. The cyclotron emission spectrum is best identified from
continuum circular polarisation. Circular spectropolarimetry allows reliable
identification of cyclotron harmonic emission and therefore determination of
the magnetic field strength. High time resolution imaging circular
polarimetry, on the other hand, can be used to map the location(s), size(s),
shape(s) and magnetic polarity of the cyclotron emission region(s) in detail.

We have carried out circular polarimetric observations at the 2.56 m Nordic
Optical Telescope (NOT) in order to detect and measure the strength of
cyclotron emission from these systems. We now proceed to describe the
observations, data analysis and results from our campaign.

Figure 1: Circular spectropolarimetry of 2MASS J0129+6715 over 1.3 orbital
cycles (left flux, right circ. polarisation). The corresponding orbital phases
calculated from the photometric ephemeris of Krushinsky et al. (2020) are
shown. The polarisation vanishes later during the first cycle than during the
second (e.g. compare orbital phases 0.533$-$0.562 vs. 1.511).

##  2 Observations

The NOT observations were mainly carried out during two observing runs in
October 2020 and March 2021 with some additional ‘fill-in’ data taken in
December 2020 (see Table 2 for the observing log). The observations were
obtained using ALFOSC, a multimode imager/spectrograph capable of both
polarimetric imaging and spectropolarimetry. Brighter systems were observed in
circular spectropolarimetric mode, whilst fainter systems were subject to
imaging circular polarimetry in single passband ($R$ band) selected to
maximize the detection probability for a typical field strength. The details
of different modes are also listed in Table 1. ALFOSC was equipped with the
e2V back illuminated deep depletion 2048x2064 CCD with 15 $\mu$m pixel size
resulting in a full FOV of 6.4x6.4’. However, only a partial window mode,
together with on-chip binning, was used in this study. The spectra were taken
using either grism #5 or #10 (see table 1.). These grisms provide spectral
coverages of 5000-10700 Å and 3300-11000 Å respectively. The grisms have
formal resolutions of $R=415$ and $R=105$ using a 1 arcsec slit. However, in
our case, we used a 1.8" polarimetric slitlet and additionally binned the CCD
by a factor of two in both directions. Thus our effective resolution is of the
order of $R\sim 50-200$ at maximum. The polarimetric observing modes involve
using a calcite block together with a rotatable $\lambda$/4 plate, rotated
with steps of 90∘. Obtaining data at four different angles of the retarder
plate allows for automatic cancelling of transmission effects in the optical
system, both for spectropolarimetic and polarimetric imaging data. However,
using only two angles (0∘ and 90∘) is also possible and provides better time
resolution for variable sources.

All the data were bias-subtracted and the spectra were flatfielded using the
halogen lamp exposures. The circular polarisation imaging time series data
were obtained under autoguiding with a quarterwave-plate rotating at 90∘
steps. Flatfielding is not required for this mode. We checked for systematic
circular polarisation by observing a bright nearby field star both in imaging
and in spectropolarimetric modes. This yielded instrumental circular
polarisation of 0.248$\pm$0.064%. We also combined 49 spectra of SDSS
J0750+4943 observed in the spectropolarimetric mode and measured the circular
polarisation within the spectral region 6500-9000Å that does not show any
cyclotron emission. We found the instrumental polarisation to be
0.052$\pm$0.016%. Compared to the observed changes of several, or even dozens
of, per cent present in our targets, the level of instrumental polarisation is
insignificant. Furthermore, during imaging polarimetry observations, the FOV
contained field stars as a check for both circular polarisation zero level and
for photometric comparison. We used IDL scripts for both the data reduction
and optimal extraction of spectra (OPTSPECEXTR package111J. Harrington,
https://physics.ucf.edu/~jh/ast/software/optspecextr-0.3.1/). Imaging aperture
photometry was carried out using our own scripts based on IDL astronomy
library.

##  3 Analysis & discussion of individual sources

In this section we proceed to describe each target and their related datasets
separately. We provide a brief introduction for each target, followed by the
description of the data as well as the relevant modelling and results. We then
give an overview of the results and discuss them jointly in the following
section.

###  3.1 2MASS J01294349+671530

2MASS J012943+6715 is a very recently discovered eclipsing 7.15 h period pre-
CV, where the secondary does not quite fill its Roche lobe (Krushinsky et al.,
2020). The eclipses of the WD show a very unusual excess in $R$ band when
compared to the other bands. It was suggested by Krushinsky et al. (2020) that
this could be explained as a result of cyclotron emission from the magnetic WD
accreting from a stellar wind with a cyclotron hump centered at/near the $R$
band.

We have obtained a set of circular spectropolarimetry data of the system
covering 8.75 h (with a 2.5 h gap during the run). The spectra cover a
wavelength range from 5500 to 8500 ÅṪhis produced 31 circular polarisation
spectra which are shown in Fig 1. Unfortunately the 2.5 h gap in the data
meant that altogether only about 0.6 of the orbital phase was covered (twice).
There is a clear cyclotron emission hump centered at around 6400-6500 Å
verifying Krushinsky et al. (2020)’s explanation for the unusual brightness of
the WD in the R-band. No other humps are visible in the 5500-10000 Å
wavelength range. This suggests that the hump is of low harmonic order. If we
identify the hump as the 2nd harmonic of the cyclotron base frequency then the
implied magnetic field strength is 81-85 MG (see eg. Wickramasinghe & Ferrario
(2000) for magnetic field determination based on cyclotron harmonics). This
would imply the base frequency cyclotron hump at around 12900 Å. However, if
we choose the 3rd harmonic instead, the magnetic field would be 53-57 MG and
the 2nd harmonic would be centered at 9740 Å. As the spectrum becomes quite
noisy towards the red end, we cannot exclude the existence of the second
harmonic in the spectra, but given these data, we prefer the 81-85 MG field
interpretation.

As our data overlaps a fraction of the orbital period (0.6), we can also
examine the level of synchronicity of the WD spin period with the binary
orbital period. The individual circularly polarised spectra show at least some
evidence for the cyclotron hump until orbital phase 0.59 during the first
orbit. However, during the second orbital cycle, the cyclotron emission
disappears between orbital phases 0.48-0.51 (Fig 1). There is thus some
evidence that the WD spin period could be $\sim$10% shorter than the orbital
period. It is likely though, that the asynchronicity would have been detected
by Krushinsky et al. (2020), even if they do not have any polarimetric data.
There are seven asynchronous polars known (Rea et al., 2017), most of which
have the WD spin and orbital periods within 1-2 %. However, there are three
systems i.e. RX J0838.7-2827 (Rea et al., 2017), Paloma (Schwarz et al., 2007;
Joshi et al., 2016) and IGR J19552+0044 (Bernardini et al., 2013), where the
WD spin period is of the order of 80-90% of the orbital period. However, none
of these have an orbital period comparable to 2MASS J0129+6715 which would
make it very unusual amongst the magnetic CVs if it too was highly
asynchronous. Further polarimetry is required to resolve this possibility.

It is interesting to note that it is not immediately clear how a system like
2MASS J0129+6715 could be produced by the evolutionary model of Schreiber et
al. (2021). While the WD appears cool enough to have started crystallization
and a brief period of asynchronous rotation is expected in this model, the
problem is that at the long period of this system, 7.15h, the donor stars in
CVs are far less over-inflated than in shorter period systems. Therefore, when
the magnetic field emerges from the WD and detaches the system, these longer
period systems will very rapidly come back into contact, since the donor star
will not substantially shrink within it’s Roche lobe. To reach the measured
Roche lobe filling factor of only 0.86 would be very difficult to achieve at
this orbital period (see supplementary figure 1 in Schreiber et al. 2021 for
example).

###  3.2 ZTF J014635.73+491443.1

Figure 2: The SED of ZTF J0146+4914 from the Vizier catalogue
(doi.org/10.26093/cds/vizier). There is no indication of any companion star.
Figure 3: Circular spectropolarimetry of ZTF J0146+4914. The left column shows
the flux spectra vs time, offset in Y-direction. The right column shows the
corresponding circular polarisation. The temporal separation of the flux
spectra is half of that of polarimetry. The labels on the right side of some
polarimetric spectra indicate the orbital phase (arbitrary 0.0 phase) coverage
and the numbers inside the panel the relevant flux spectra they are based on.
The orbital modulation of cyclotron humps is clearly visible both in the flux
spectra and the circular polarisation spectra. Figure 4: Circular $R$ band
imaging photopolarimetry of SDSS J0308-0054. Top panel: the differential
magnitude compared to a nearby comparison star. Lower panel: the circular
polarised light curve.

ZTF J0146+4914 was discovered as a candidate polar with a photometric period
of 2.057 h by Guidry et al. (2021). The system shows cyclotron hump features
in the optical spectra, compatible with a WD magnetic field of $\sim$89 MG.
The distance of (56.3 $\pm$ 0.3 pc, Bailer-Jones et al. (2018)) could
potentially establish ZTF J0146+4914 as the closest known polar. However,
there are no clear emission lines visible in the spectra, suggesting a wind
accreting system. Furthermore, there is also no clear sign of the secondary in
the SED Fig 2, indicating a very low mass companion, probably a brown dwarf.
This, together with the orbital period, makes the system a very likely period
bouncer (i.e. a system that has passed the period minimum in its orbital
evolution and is now evolving towards longer periods) and thus a very old
system.

We have attempted to characterise the donor brown dwarf by extracting the WISE
(Wright et al., 2010) magnitudes from the CatWISE (Eisenhardt et al., 2020)
catalogue and comparing them against the theoretical DA WD WISE magnitudes
obtained by interpolating from the Montreal WD model grid (Holberg & Bergeron,
2006; Tremblay, Bergeron, & Gianninas, 2011)
222http://www.astro.umontreal.ca/~bergeron/CoolingModels using $T=8700$K and
$\log(g)=8.355$ from Gentile Fusillo et al. (2021). We detect a very small IR
excess in the WISE bands W1 and W2. The corresponding absolute WISE W2
magnitude of 13.65 would imply a $\sim$T8 donor based on Fig. 29 of
Kirkpatrick et al. (2011). However, the colour of the IR excess, $W1-W2=0.96$,
suggests a much earlier donor of T2-T3 spectral class (Fig. 1 of Kirkpatrick
et al. 2011). It is possible that the system contains a late T type donor,
that has been slightly heated up by its WD companion to match these figures.

In fact, ZTF J0146+4914 appears to belong to a distinct subpopulation of
systems, which include so far EF Eri (Williams et al., 1979), SDSS J1212+0136
(Schmidt et al., 2005b), IL Leo (Schmidt et al., 2007; Parsons et al., 2021)
and the SDSS J1250+1549 & SDSS J1514+0744 (Breedt et al., 2012). All of these
share the following properties: their donor stars are very likely brown
dwarfs, they show prolonged states of low accretion and they have orbital
periods below the period gap ($\la$2hrs). Furthermore, all of them also
exhibit WD temperatures of the order of 10 000 K (ZTF J0146+4914 appears in
the eDR3 WD catalogue of Gentile Fusillo et al. 2021 with a WD effective
temperature of 8700K). The magnetic field strengths of these systems vary from
13 up to 89 MG; effectively covering the full range typically observed in
polars. These systems may represent highly evolved polars, i.e. the equivalent
of WZ Sge dwarf novae among the non-magnetic CVs. The fact that one of these
system (ZTF J0146+4914) is located only 56 pc away suggests that these systems
are likely to be fairly common, as was already suggested by Patterson (2011)
and Breedt et al. (2012) based on their space density estimates.

The previously detected cyclotron humps in the optical spectra of ZTF
J0146+4914 are centered at $\sim$4000 Å and $\sim$6000 Å and show somewhat
different behaviour as a function of orbital phase. Especially the hump at
6000 Å, assumed to be the 2nd harmonic of the cyclotron base frequency,
changes remarkably in both width and shape over the orbital period.

We obtained circular spectropolarimetry of ZTF J0146+4914 at two epochs in
December 2020 using ALFOSC with grism #10 covering a wavelength range of
3500-8500Å. We show the higher time resolution spectropolarimetric time series
in Fig 3. There are a couple of striking features in the spectra. Firstly, we
detect the changes in width and shape of the 6000 Å cyclotron hump over the
orbital period, discovered by Guidry et al. (2021). Secondly, we note that the
sign of the circular polarisation is reversed both at the red and blue ends of
the cyclotron hump. This reversal persists for the duration the cyclotron
region is visible over the orbital period, suggesting that there are in fact
two emission regions visible at the same longitude on the WD surface. These
would have different polarity, but similar visibility over the orbital period.
It also appears that the cyclotron hump produced by the positive region is
wider than the one produced by the negative region, leading to the observed
sign reversal effect. Since the humps are centered at the same wavelength, the
magnetic field strengths at the two accretion regions are very similar. The
different hump widths are thus likely to originate from differences in optical
depth due to plasma density/temperatures together with the different viewing
angles (however, see our modelling attempts on source SDSS J0750+4943 in §3.4,
which shows very similar behaviour). We are not aware that such wavelength
dependent cyclotron hump sign reversal has been reported before in any
magnetic CV. For instance, there is no suggestion of this in the
spectropolarimetry of LARPs in Schmidt et al. (2005a) or Schmidt et al.
(2007).

###  3.3 SDSS J030856.55-005450.7

This source is an eclipsing dM4+WD binary with a 4.5 h period (Becker et al.,
2011) . There is some tentative evidence that the eclipse depth may be
variable (from CCD photometry with 0.5m class telescopes), which might
indicate the presence of changing accretion rate/cyclotron emission on the WD
surface. We first attempted circular spectropolarimetry, but the source was
too faint to obtain useful data. The source was then observed for 4.2 h in
imaging circular polarisation mode (Fig 4) under variable transparency
conditions. The $R$ band photometry reveals the orbital modulation, but there
is no indication of any circular polarisation.

###  3.4 SDSS J075015.11+494333.2

Figure 5: Circular spectropolarimetry of SDSS J0750+4943. The four top panels
show the flux and circular polarisation of spectrum #1 (left) and #8 (right)
highlighting the orbital phases where the cyclotron emission regions are
visible and hidden behind the WD. The bottom four panels show the trailed
cyclotron emission spectra and its circular polarisation (top two), together
with our cyclotron model fits (bottom two). Note that there is a 0.5 orbital
period gap between spectra # 7 and 11. The flux and polarisation within this
gap also mark the zero levels accordingly.

SDSS J075015.11+494333.2 has an orbital period of 4.2 h (Drake, 2014) and was
identified as a pre-polar candidate from its unusual light curve (Parsons et
al., 2015). Optical spectra revealed strong cyclotron features making the
source a very likely high field polar, with a magnetic field of either 99 or
196 MG (Parsons et al., 2021). We obtained 4.4 h of circular
spectropolarimetry using a grism #5 with ALFOSC on 4th Mar 2021. The
polarimetric time resolution was 600 s resulting in 26 spectra over the
orbital period (Fig. 5). The spectra cover a wavelength range of 5000-9500Å.
We detect a single cyclotron hump in this range, centered at around
5600-5700Å. There is no indication of circular polarisation towards the red
end of the spectrum. This implies that the detected cyclotron hump must
correspond to the 2nd harmonic of the cyclotron base frequency with a magnetic
field of 94-96 MG, consistent with the VLT/X-shooter results of Parsons et al.
(2021).

The circular polarisation behaviour of the cyclotron hump is intriguing in a
similar way as in ZTF J0146+4916 i.e. the polarity shows sign reversal at both
blue and red ends of the cyclotron hump. Again, like in ZTF J0146+4914, both
the positive and negative poles seem to show the same orbital behaviour. Thus
the same conclusions regarding the accretion geometry/physics apply here. As
the observations of SDSS J0750+4943 are of much higher quality than the ZTF
J0146+4914 data, we have attempted modelling of the spectropolarimetry using a
simple model, consisting of two cyclotron emission regions on the WD surface.
Both of these are assumed to have a “slab" like geometry with radial magnetic
field lines. The cyclotron emission spectrum, dictated by the viewing angle,
surface magnetic field strength, the constant plasma temperature and a
constant dimensionless plasma parameter $\Lambda$, is computed using the
method (and code) described in Woelk & Beuermann (1992),Woelk & Beuermann
(1993) and Rousseau et al. (1996).

In order to fit the spectropolarimetry, we have first taken advantage of the
fact that the cyclotron emission is not visible at all orbital phases. We can
thus subtract away the contribution from the companion star from all the
spectra showing evidence for cyclotron emission. We have then fitted
simultaneously 24 flux spectra and the corresponding circular polarisation
spectra from different orbital phases. The model includes two point-like
cyclotron regions that are free to move on the WD surface. In addition to
their locations, the free parameters for the two cyclotron regions
(separately) include: The magnetic field strength, shock temperature,
dimensionless plasma parameter $\Lambda$ and relative brightness of the
regions. We attempted fitting with both: i) different values of fixed
inclination and ii) the inclination as an additional free parameter. The
fitting itself was carried out using a Differential Evolution (DE) algorithm
(Storn & Price, 1997). The results are shown in Fig 5. Whilst we are able to
obtain reasonable fits to the flux spectrum, the fits to the circular
polarisation spectrum are not good. Furthermore, the resulting models place
both of the cyclotron emission regions in the ‘southern’ hemisphere of the WD
(colatitudes depending on the inclination), separated only by a couple of
degrees. Yet these two regions are required to have different magnetic
polarity. In fact, the flux spectrum itself could be reasonably well fitted
with just a single simple cyclotron region, if the polarisation would be
ignored. Clearly, our simple model is not likely to represent the true picture
and more detailed modelling would be required, perhaps employing a multipole
magnetic field for the WD.

###  3.5 SDSS J085336.03+072033.5

SDSS J0853+0720 shows a 3.6 h period (Nebot, 2011), with X-Shooter spectra
revealing cyclotron features suggesting a magnetic WD with a field strength of
84 MG (Parsons et al., 2021). Our $R$ band photopolarimetric time series spans
over a single 4 h run (Fig 6). The light curve shows a clear overlap with
repeating behaviour in the end. It is characterised by two maxima of unequal
brightness. The lower maximum coincides with a clear maximum in negative
circular polarisation, which reaches about $\sim$8%. The circular polarisation
is reversed during the higher maximum but, somewhat unexpectedly, the circular
polarisation only reaches about +3 % and lasts for much less than the duration
of the hump in the light curve. Nevertheless, we suggest this system is
accreting on two poles, the visibility of which alternates over the orbital
period. Finally, the $R$ band’s effective wavelength of 6400 Å coincides with
the centre wavelength of the second harmonic cyclotron hump from the 84 MG
field.

Figure 6: Circular $R$ band photopolarimetry of SDSS J0853+0720. Top panel:
the differential magnitude compared to a nearby comparison star. Lower panel:
the circular polarised light curve. Figure 7: Circular photopolarimetry of
SDSS J1212+0136. Top panel: the differential magnitude compared to a nearby
comparison star. Lower panel: the circular polarised light curve. Figure 8:
Circular spectropolarimetry of SDSS J2229+1853, The panels (from top) show the
count spectrum, the circular polarisation spectrum, the polarised counts and
the unpolarised count spectrum, associated mostly with the M dwarf. A very
prominent cyclotron hump centered at around 5800-6000 Å is present.

Figure 9: Top panels: Circular photopolarimetry of SDSS J2229+1853 with a
cyclotron emission model fit. Lower panel: the implied cyclotron emission map
of the WD surface. The viewing angle to the WD is shown as at $i$=45,
phase=0.0

###  3.6 SDSS J121209.31+013627.7

This is a very likely brown dwarf + magnetic WD binary with an orbital and WD
spin period of 88 min (Schmidt et al., 2003, 2005b; Burleigh et al., 2006).
Stelzer et al. (2017) confirmed it as weakly accreting polar-like object
through time-resolved X-ray observations with XMM-Newton. A WD magnetic field
of 7 MG (mean) and 13 MG (dipole) were reported by Schmidt et al. (2003) and
supported based on NIR spectroscopy (Farihi, Burleigh, & Hoard, 2008). Given
the low magnetic field estimate, we obtained $i$ band data of this system
(Fig.7). However, the observations were cut short due to bad weather and the
data cover only 52 min (i.e. 0.59 in orbital phase). We do not see any
evidence for circular polarisation, but it is possible the polarisation signal
was missed due to the limited phase coverage. Furthermore, the observations
suffered from thin cirrus.

###  3.7 SDSS J222918.95+185340.2

SDSS J222918.95+185340.2 was identified as a potential pre-polar candidate
from its unusual CRTS light curve (Drake, 2014). X-shooter spectra covering
2/3 of the orbital cycle reveal at least 6 cyclotron lines appearing at
different phases, suggesting a complicated field structure (Parsons et al.,
2021). Optical light is very dominated by the M3 companion star in the very
red end of the spectrum. The binary star’s orbital period is 4.5 h and the
system is extremely close to Roche-lobe filling (the fit is actually
consistent with Roche-lobe filling) but beyond the cyclotron lines there are
no other obvious signs of accretion (Parsons et al., 2021).

We obtained 4 circularly polarised spectra of SDSS J2229. The summed up flux
and circular polarisation spectra are shown in Fig 8. We detect a clear
cyclotron hump centered around 6000 Å. There is no evidence in circular
polarisation for any other cyclotron hump up to 10000 Å. This strongly
suggests that the 6000 Å hump is likely the second harmonic component, arising
from the 85-90 MG WD magnetic field, consistent with the X-shooter data
(Parsons et al., 2021).

The $R$ band circular time resolved photopolarimetry is shown in Fig 9. The
circular polarisation reaches 30% at its peak, compatible with the
spectropolarimetric results (Fig 8). Perhaps the most striking feature of the
photopolarimetry is the total absence of any flickering in the data. This is
likely due to wind accretion, as there is no evidence for any accretion stream
or Roche lobe overflow in the system (i.e. no emission lines or signs for
accretion stream continuum emission). It is feasible that accreted stellar
wind lacks the clumpiness of accretion flow resulting from Roche lobe
overflow.

We have used the resulting high signal-to-noise flux and polarisation light
curves to map the WD surface in cyclotron radiation. This approach is very
similar to the ’Stokes Imaging’ (Potter, Hakala, & Cropper, 1998), although we
are lacking the linear polarisation information. The approach is also related
to the one used for TESS CD Ind data (Hakala et al., 2019), although the TESS
data obviously lacked any polarisation information.

In brief, we have utilised the cyclotron emission models of Wickramasinghe &
Meggitt (1985) to model the dependence of the cyclotron flux and circular
polarisation as a function of angle between the line of sight and the magnetic
field line. We have used a 10 keV, 2nd harmonic model with $\Lambda=10^{5}$
for this. Since we do not have information on the magnetic field geometry of
the WD, we have assumed that the cyclotron emission originates near the
magnetic poles, where the magnetic field lines approximately align with the
surface normal. We have also assumed that the cyclotron emission is not
confined to the surface of the WD, but the emission regions can have height up
to 6% of the WD radius (corresponding to observing the shock front up to 20∘
beyond the limb crossing). This is directly indicated by the data, as circular
polarisation sign reversals are seen before and after the limb crossings (i.e.
the cyclotron emission region is seen from below). As there is no evidence for
an obvious second pole (with negative circular polarisation), we have assumed
that all the field lines, where the emission is detected, flow radially
outwards. However, we have also experimented with sign reversal of magnetic
field for different hemispheres, but that failed to produce acceptable fits.
For the actual fitting, we used a Differential Evolution (DE, Storn & Price
1997) algorithm to optimize a smoothest possible cyclotron emission map over
the WD surface. This was done using a grid with 4∘ resolution (see Hakala et
al. 2019 for more details). The fitting was carried out with fixed inclination
angle. Several inclination angles were tried, together with varying amount of
unpolarised additional flux.

The best results were obtained using $i=45^{\circ}$ and no additional
unpolarised light. This strongly suggests there is no contribution from an
accretion stream which, if present, typically accounts for large fraction of
optical emission, as evidenced by eclipse profiles (Hakala, 1995) and Doppler
mapping (Heerlein, Horne, & Schwope, 1999). The resulting cyclotron map is
shown in the bottom panel of Fig 9. The cyclotron emission appears to be
dominated by a main accretion region in the ‘Northern’ hemisphere at a
latitude of 40-50∘ and pointing almost towards the secondary star. There are
two additional weaker regions roughly at the same WD longitude, but at
latitude -20∘. These secondary regions are thus separated by $\sim$60-70∘ from
the main region. Given the small gap between these regions, the separation
into two close regions is probably not significant, but likely arises from the
regularisation. Remarkably, the locations of the main cyclotron regions Fig 9
(at different WD hemispheres, both facing the donor star) are almost identical
to the ones measured in another similar system SDSS J030308.35+005444.1
(Parsons et al., 2013) using eclipse data, even if the magnetic field here is
an order of magnitude greater. However, we should be careful in putting too
much weight on interpreting these results though, as the true magnetic field
configuration is not known nor implemented in our modelling, which assumes
radial field lines at cyclotron emission sites.

Source | POrb | Details | Peak circ.pol. | Inferred/known B-field | Remarks
---|---|---|---|---|---
2MASS J01294349+6715301 | 7.15 h | K7V+WD, ecl. | 4% | 81-85 MG (or 53-57 MG) | 53-57 MG field cannot be ruled out.
ZTF J014635.73+491443.14 | 2.057 h | BD?+WD |  $\sim$40% | 89 MG4 | Distance 56.3 $\pm$ 0.3 pc
SDSS J030856.55$-$005450.73 | 4.5 h | M4+WD, ecl. | No detection | N/A |
SDSS J075015.11+494333.22 | 4.2 h | M2.5+WD |  $\sim$20% | 94-96 MG2 | Complex magnetic field geometry implied.
SDSS J085336.03+072033.52 | 3.6 h | M4+WD | 8% | 84 MG2 | two pole accretion.
SDSS J121209.31+013627.75 | 1.47 h | BD+WD | No detection | 13 MG5 | Observations covered only 0.59 in phase
SDSS J222918.95+185340.22 | 4.5 h | M3+WD | 45% | 85-90 MG | Cyclotron mapping reveals two-pole accretion.
Table 2: The summary of system parameters and main results. .(1Krushinsky et
al. (2020),2Parsons et al. (2021),3Becker et al. (2011),4Guidry et al.
(2021),5Schmidt et al. (2003) Figure 10: The cyclotron emission spectra from
constant low temperature (1keV) and ($log\Lambda=5$) “slab" models as a
function of the magnetic field strength. The locations of g and r passbands
are shown with vertical dashed lines.

##  4 Conclusions

We have obtained the first optical circular imaging and spectropolarimetry of
seven pre-polar (or low accretion rate polar) candidates. The sources were
chosen based on either existing optical spectroscopy or optical light curve
shapes suggesting cyclotron emission in the systems. None of the chosen
systems show strong emission lines or any other evidence of accretion taking
place and light curve modelling of three of these systems (Parsons et al.,
2021) indicate that the secondary stars do not fill their Roche lobes. We find
that in 5 out of 7 systems we detect variable optical circular polarisation,
confirming the cyclotron emission and thus the presence of the magnetic WD.
The summary of our polarimetric results can be found in Table 2.

Intriguingly, all of the 5 systems exhibiting circular polarisation show WD
magnetic fields in excess of 80 MG. This is considerably higher than is
typically found in Roche lobe overflow accreting polars (Ferrario, de Martino,
& Gänsicke, 2015). In fact Ferrario, de Martino, & Gänsicke (2015) only list
four polars (or pre-polars) out of 72 that show WD field strength in excess of
80 MG. This suggests that maybe the pre-polar stage systems with the highest
field WDs are more easily detected than the lower field pre-polars. This could
be due to the more effective channeling of the stellar wind from the secondary
star on to the WD. The second option is that the WDs in pre-polars indeed are
intrinsically divergent and have experienced a different formation and/or
accretion history than the WDs in polars. Finally, the most likely case is
perhaps that this is due to a selection effect. Schwope et al. (2009) discuss
the similar difference in WD properties regarding their sample or 9 LARPs that
showed higher then average magnetic fields, cool WDs and longer orbital
periods than the known polars. They conclude that selection effects are to
likely play a significant role, as the accreting plasma in pre-polars is
tenuous and cooler than in polars, that sustain much higher accretion rates.
This leads to the cyclotron emission peaking in the lower harmonics and only
weak emission from harmonics higher than 3rd is observed.

We demonstrate the effect in Fig 10, where we have plotted the cyclotron
spectra from constant temperature models (log $\Lambda$=5 and Teff = 1 keV,
viewing angle = 45∘) with magnetic field strength ranging from 40 to 100 MG.
The two strong cyclotron humps correspond to the 2nd and 3rd harmonics. We
also show the Sloan g and r passbands separated by vertical dashed lines. It
is clear that the systems containing more magnetic WDs are much more easily
picked up based on say Sloan passband photometry and/or optical spectroscopy.
The same bias is not equally obvious for the polars, where the shock region
near the WD surface is much hotter and denser, leading to stronger cyclotron
emission at harmonics higher than the 3rd. Thus it is likely that the range of
WD magnetic fields in our sample of systems is at least partly due to the fact
that our sample is selected based on possible cyclotron features in the
spectra or light curves of the candidates.

It is worth mentioning that, even if the fundamental frequency of cyclotron
emission is generally thought to be optically thick and non-polarised in
polars, this might not be the case in wind-accreting systems, where the
accretion rates are several magnitudes lower Schwope et al. (2002). Thus it is
feasible, based on our cyclotron model computations, that the emission at the
fundamental frequency becomes optically thin and exhibits significant circular
polarisation if the plasma density becomes low enough (log$\Lambda\lesssim$
2.5 in our modelling). The implication is that the magnetic field estimates
based on seeing a single cyclotron hump, with no polarisation towards the red
end of the optical spectrum, could be too low by a factor of two, in case the
observed hump would in fact represent emission at the fundamental frequency.
This applies to sources 2MASS J0129+6715, SDSS J0750+4943 and SDSS J2229+1853.

We have attempted to model the cyclotron emission in detail in two of the
systems. The spectropolarimetric modelling of SDSS J0750+4943 failed to
produce satisfactory results suggesting a complex WD field geometry. However,
the photopolarimetric modelling of SDSS J2229+1853 was successful and implied
two elongated cyclotron emission region located roughly at the same
longitudes, but with 70-80∘ separation in latitude. The light curve is
entirely dominated by the cyclotron emission, with no additional emission
whatsoever. It is worth noting that Parsons et al. (2021) measured the
accretion rate for SDSS J2229+1853 (based on cyclotron flux). They found
$6.5\times 10^{-13}$M⊙ yr-1, which provides the lower limit for the stellar
wind from the M3V companion star. This value is very high
($30\dot{M}_{\odot}$), but such winds have been measured (Wood et al., 2021).
However, the likely implication is that almost all of the stellar wind is
accreted and, at least in its current state, the magnetic braking should be
very inefficient. The total absence of any flickering on top of the cyclotron
emission is a strong argument against the typical accretion scheme in polars
and perhaps serves as an indication of accretion via stellar wind.

Apart from 2MASS 012943+6715, and ZTF J0146+4914, all the other three systems
containing a highly magnetic WD are compatible with the pre-polar scheme as
outlined in Schreiber et al. (2021). The suggested asynchronism in 2MASS
012943+6715 needs to be verified by further polarimetric observations.
Finally, we detect a very likely brown dwarf companion in the period bouncer
system ZTF J0146+4914. We estimate that the companion brown dwarf is possibly
of late T spectral type, slightly heated by the WD, to match the WISE
observations.

## Acknowledgements

Based on observations made with the Nordic Optical Telescope, owned in
collaboration by the University of Turku and Aarhus University, and operated
jointly by Aarhus University, the University of Turku and the University of
Oslo, representing Denmark, Finland and Norway, the University of Iceland and
Stockholm University at the Observatorio del Roque de los Muchachos, La Palma,
Spain, of the Instituto de Astrofisica de Canarias. The data presented here
were obtained with ALFOSC, which is provided by the Instituto de Astrofisica
de Andalucia (IAA) under a joint agreement with the University of Copenhagen
and NOT. We would like to thank the NOT staff for carrying out the
observations. Armagh Observatory and Planetarium is core funded by the
Northern Ireland Executive through the Dept. for Communities. SGP acknowledges
the support of a Science and Technology Facilities Council (STFC) Ernest
Rutherford Fellowship.213056.71+442046.5. TRM and BTG were supported by grant
ST/T000406/1 from the Science and Technology Facilities Council (STFC). This
project has received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (Grant
agreement No. 101020057). This research has made use of the VizieR catalogue
access tool, CDS, Strasbourg, France (DOI : 10.26093/cds/vizier). The original
description of the VizieR service was published in 2000, A&AS 143, 23. For the
purpose of open access, the author has applied a Creative Commons Attribution
(CC BY) licence to any Author Accepted Manuscript version arising. Finally, we
wish to thank the referee, Prof. John Landstreet, for detailed remarks that
certainly improved the paper considerably.

## Data Availability

The data will be shared on reasonable request to the corresponding author.

## References

  * Bailer-Jones et al. (2018) Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Mantelet G., Andrae R., 2018, AJ, 156, 58
  * Becker et al. (2011) Becker A. C., Bochanski J. J., Hawley S. L., Ivezić Ž., Kowalski A. F., Sesar B., West A. A., 2011, ApJ, 731, 17.
  * Belloni et al. (2020) Belloni D., Schreiber M. R., Pala A. F., Gänsicke B. T., Zorotovic M., Rodrigues C. V., 2020, MNRAS, 491, 5717.
  * Bernardini et al. (2013) Bernardini F., de Martino D., Mukai K., Falanga M., Andruchow I., Bonnet-Bidaud J.-M., Masetti N., et al., 2013, MNRAS, 435, 2822.
  * Breedt et al. (2012) Breedt E., Gänsicke B. T., Girven J., Drake A. J., Copperwheat C. M., Parsons S. G., Marsh T. R., 2012, MNRAS, 423, 1437. doi:10.1111/j.1365-2966.2012.20965.x
  * Burleigh et al. (2006) Burleigh M. R., Marsh T. R., Gänsicke B. T., Goad M. R., Dhillon V. S., Littlefair S. P., Wells M., et al., 2006, MNRAS, 373, 1416.
  * Cropper (1990) Cropper M., 1990, SSRv, 54, 195.
  * Drake (2014) Drake, A. J., et al., 2014, ApJS, 213, 9
  * Eisenhardt et al. (2020) Eisenhardt P. R. M., Marocco F., Fowler J. W., Meisner A. M., Kirkpatrick J. D., Garcia N., Jarrett T. H., et al., 2020, ApJS, 247, 69. doi:10.3847/1538-4365/ab7f2a
  * Farihi, Burleigh, & Hoard (2008) Farihi J., Burleigh M. R., Hoard D. W., 2008, ApJ, 674, 421.
  * Ferrario, de Martino, & Gänsicke (2015) Ferrario L., de Martino D., Gänsicke B. T., 2015, SSRv, 191, 111.
  * Gentile Fusillo et al. (2021) Gentile Fusillo N. P., Tremblay P.-E., Cukanovaite E., Vorontseva A., Lallement R., Hollands M., Gänsicke B. T., et al., 2021, MNRAS, 508, 3877. doi:10.1093/mnras/stab2672
  * Guidry et al. (2021) Guidry J. A., Vanderbosch Z. P., Hermes J. J., Barlow B. N., Lopez I. D., Boudreaux T. M., Corcoran K. A., et al., 2021, ApJ, 912, 125.
  * Hakala (1995) Hakala P. J., 1995, A&A, 296, 164.
  * Hakala et al. (2019) Hakala P., Ramsay G., Potter S. B., Beardmore A., Buckley D. A. H., Wynn G., 2019, MNRAS, 486, 2549.
  * Heerlein, Horne, & Schwope (1999) Heerlein C., Horne K., Schwope A. D., 1999, MNRAS, 304, 145
  * Holberg & Bergeron (2006) Holberg J. B., Bergeron P., 2006, AJ, 132, 1221. doi:10.1086/505938
  * Isern et al. (2017) Isern J., García-Berro E., Külebi B., Lorén-Aguilar P., 2017, ApJL, 836, L28. doi:10.3847/2041-8213/aa5eae
  * Joshi et al. (2016) Joshi A., Pandey J. C., Singh K. P., Agrawal P. C., 2016, ApJ, 830, 56.
  * Kirkpatrick et al. (2011) Kirkpatrick J. D., Cushing M. C., Gelino C. R., Griffith R. L., Skrutskie M. F., Marsh K. A., Wright E. L., et al., 2011, ApJS, 197, 19. doi:10.1088/0067-0049/197/2/19
  * Krushinsky et al. (2020) Krushinsky V., Benni P., Burdanov A., Antokhin I., Antokhina E., Jehin E., Barkaoui K., et al., 2020, MNRAS, 493, 5208.
  * Li, Wu, & Wickramasinghe (1994) Li J. K., Wu K. W., Wickramasinghe D. T., 1994, MNRAS, 268, 61.
  * Nebot (2011) Nebot Gómez-Morán A., et al., 2011, A&A, 536, A43
  * Paczyński (1967) Paczyński B., 1967, AcA, 17, 287 Paczyński, 1976, in Eggleton P, Mitton S, Whelan J (eds) Structure and evolution of close binary systems. IAU symposium, vol 73, p 75.
  * Parsons et al. (2013) Parsons S. G., Marsh T. R., Gänsicke B. T., Schreiber M. R., Bours M. C. P., Dhillon V. S., Littlefair S. P., 2013, MNRAS, 436, 241.
  * Parsons et al. (2015) Parsons S. G., et al., 2015, MNRAS, 449, 2194
  * Parsons et al. (2021) Parsons S. G., Gänsicke B. T., Schreiber M. R., Marsh T. R., Ashley R. P., Breedt E., Littlefair S. P., et al., 2021, MNRAS, 502, 4305.
  * Patterson (2011) Patterson J., 2011, MNRAS, 411, 2695. doi:10.1111/j.1365-2966.2010.17881.x
  * Potter, Hakala, & Cropper (1998) Potter S. B., Hakala P. J., Cropper M., 1998, MNRAS, 297, 1261.
  * Rea et al. (2017) Rea N., Coti Zelati F., Esposito P., D’Avanzo P., de Martino D., Israel G. L., Torres D. F., et al., 2017, MNRAS, 471, 2902.
  * Rousseau et al. (1996) Rousseau T., Fischer A., Beuermann K., Woelk U., 1996, A&A, 310, 526
  * Schmidt et al. (2003) Schmidt G. D., Harris H. C., Liebert J., Eisenstein D. J., Anderson S. F., Brinkmann J., Hall P. B., et al., 2003, ApJ, 595, 1101. doi:10.1086/377476
  * Schmidt et al. (2005a) Schmidt G. D., Szkody P., Vanlandingham K. M., Anderson S. F., Barentine J. C., Brewington H. J., Hall P. B., et al., 2005, ApJ, 630, 1037. doi:10.1086/431969
  * Schmidt et al. (2005b) Schmidt G. D., Szkody P., Silvestri N. M., Cushing M. C., Liebert J., Smith P. S., 2005, ApJL, 630, L173. doi:10.1086/491702
  * Schmidt et al. (2007) Schmidt G. D., Szkody P., Henden A., Anderson S. F., Lamb D. Q., Margon B., Schneider D. P., 2007, ApJ, 654, 521. doi:10.1086/509613
  * Schreiber et al. (2021) Schreiber M. R., Belloni D., Gänsicke B. T., Parsons S. G., Zorotovic M., 2021, NatAs, 5, 648. doi:10.1038/s41550-021-01346-8
  * Schwarz et al. (2007) Schwarz R., Schwope A. D., Staude A., Rau, A., Hasinger G., Urrutia T., Motch, C., 2007, A&A, 473(2), 511. doi:10.1051/0004-6361:20077684
  * Schwope et al. (2002) Schwope A. D., Brunner H., Hambaryan V., Schwarz R., 2002, ASPC, 261, 102
  * Schwope et al. (2009) Schwope A. D., Nebot Gomez-Moran A., Schreiber M. R., Gänsicke B. T., 2009, A&A, 500, 867. doi:10.1051/0004-6361/200911699
  * Stelzer et al. (2017) Stelzer B., de Martino D., Casewell S. L., Wynn G. A., Roy M., 2017, A&A, 598, L6. doi:10.1051/0004-6361/201630038
  * Storn & Price (1997) Storn R. & Price K., 1997, Journal of Global Optimization, 11, 4, 341.
  * Tremblay, Bergeron, & Gianninas (2011) Tremblay P.-E., Bergeron P., Gianninas A., 2011, ApJ, 730, 128. doi:10.1088/0004-637X/730/2/128
  * Verbunt & Zwaan (1981) Verbunt F., Zwaan C., 1981, A&A, 100, L7
  * Wickramasinghe & Meggitt (1985) Wickramasinghe D. T., Meggitt S. M. A., 1985, MNRAS, 214, 605.
  * Wickramasinghe & Ferrario (2000) Wickramasinghe D. T., Ferrario L., 2000, PASP, 112, 873.
  * Williams et al. (1979) Williams G., Johns M., Price C., Hiltner A., Boley F., Maker S., Mook D., 1979, Natur, 281, 48. doi:10.1038/281048a0
  * Woelk & Beuermann (1992) Woelk U., Beuermann K., 1992, A&A, 256, 498
  * Woelk & Beuermann (1993) Woelk U., Beuermann K., 1993, A&A, 280, 169
  * Wood et al. (2021) Wood B. E., Müller H.-R., Redfield S., Konow F., Vannier H., Linsky J. L., Youngblood A., et al., 2021, ApJ, 915, 37.
  * Wright et al. (2010) Wright E. L., Eisenhardt P. R. M., Mainzer A. K., Ressler M. E., Cutri R. M., Jarrett T., Kirkpatrick J. D., et al., 2010, AJ, 140, 1868. doi:10.1088/0004-6256/140/6/1868