\section{Introduction} \label{s:intro}

The study of magnetism in stars is important to understand its role in the formation and evolution of stars in general (fragmentation of protostellar clouds, angular momentum lost, convection, winds, magnetic activity cycles, etc.).
Since characterizing main-sequence magnetic properties can be rather complicated due to their low field strengths, studying their remnants, mostly white dwarf stars, can provide valuable information about the magnetic characteristics of their progenitors as a whole.
Magnetic white dwarfs (MWDs) that are not in an interacting binary system represent around 13 \% of all white dwarfs.
Magnetism in these former stellar cores have been observed over a wide range of strength varying from $10^3$ to $10^9$ G.
The exact origin of the magnetic field in most white dwarf stars is still the subject of debate \citep{Ferrario2015a, Ferrario2015b, Ferrario2020}: are the fields fossil or the results of some binary interaction, merging events or internally produces as core crystallization begins?
Are there signatures in the observed properties of a MWD that allow us to infer its evolutionary channel?

In the fossil field scenario \citep{Tout2004, Ferrario2015b}, the high magnetic fields observed in white dwarf stars are simply the results of magnetic flux conservation as the star evolves (collapses) to a much smaller radius.
Magnetic main-sequence Ap and Bp stars thus appear as natural candidate as progenitors of magnetic white dwarfs (MWDs).
The correspondence between the expected vs. observed field intensities from this population is a strong argument in favor of this hypothesis.
However, the required birth rate of those main-sequence stars is insufficient by a factor of 2 $-$ 3, to account for the incidence of magnetism in WDs \citep{Wickramasinghe2005}.
Alternatively, high magnetic field could be the results of some sort of binary interaction.
The lack of binary systems involving a MWD and a fully detached (and non-degenerate) companion supports this idea \citep{Liebert2015}, at least for a fraction of the population of MWDs \citep{Tout2008}.
Interestingly, MWDs are found to have an average mass that is much higher than their non-magnetic counterparts \citep[average mass of \Msol{\sim 0.8} vs \Msol{0.6},][]{Ferrario2015a}, suggesting the merger evolution branch to be more likely.
It must be noted, however, that \citet{Bagnulo2021} recently showed, based on the local sample (40 pc) that the higher than average mass for MWDs could probably be a selection effect.
Indeed, according to their findings, high field MWDs of normal mass only appear much later on the cooling sequence (possibly as the result of a dynamo mechanism that is triggered by the start of core crystallization) while the high mass ones are present at all temperatures.
Consequently, the overluminous hot MWDs could be over-represented in magnitude limited surveys, skewing the mass distribution toward higher value.

Studying the incidence of magnetism among white dwarf stars, their distribution in terms of mass, effective temperature, field intensity and geometry should provide valuable hints to help us understand the evolutionary paths leading to MWDs.

The incidence of magnetism among white dwarf stars is still uncertain as volume-limited samples and magnitude-limited sample yield values ranging from 10 to 20\% \citep{Kawka2007, Giammichele2012, Bagnulo2021} to 2-5\% \citep{Liebert2003, Kepler2015} respectively.
However, those numbers certainly suffer, at the lowest field, from the fact that spectropolarimetric measurements have been obtained only for a small fraction of stars \citep[and references therein]{Vennes2018, Kawka2007, Kawka2019, Landstreet2012, Landstreet2015, Bagnulo2022}.
Moreover, the low signal-to-noise spectra for many objects barely allow the line splitting to be resolved \citep[and references therein]{Kepler2015, Kepler2019}.
It is to be noted that the number of known magnetic white dwarfs has increased substantially in the last two decades or so.
Indeed, fewer than 70 objects were known at the turn of the century \citep{Wickramasinghe2000} but thanks mainly to the Sloan Digital Sky Survey \citep[SDSS,][]{York2000} there are allegedly now more than 800 objects that have tentatively been identified as MWD \citep{Gansicke2002, Schmidt2003, Vanlandingham2005, Kulebi2009, Kepler2013, Kepler2015, Ferrario2015a}.
About 75\% of these are isolated stars with estimated magnetic field strength between $10^3$ G to $10^3$ MG, the other 25\% being in interacting binary systems \citep{Ferrario2020}.

However, our knowledge of the bulk properties of MWDs with magnetic fields in the MG regime relies mostly on pioneering work done in the 1980’s and 1990’s \citep[for example][]{Wickramasinghe1979, Jordan1991, Jordan1992, Bergeron1992, Wickramasinghe2000, Schmidt2003}.
The last comprehensive homogeneous analysis of a large hydrogen-rich MWD sample (100+ object), that of \citet{Kulebi2009}, is now more than a decade old.
It must be noted that most studies, including that of \citet{Kulebi2009}, have been carried out assuming $\log{g} = 8$ as no parallax measurements were available for most objects at the time of analysis.
Another reason is that, as mentionned in \citet{Ferrario2020}, a common method for the determination of the surface gravity for the non-magnetic WDs is performed using the broadening of spectral lines due to the Stark effect.
Since no calculations of line profiles taking both the Stark and magnetic broadening into account exists, determination of $\log{g}$ using the spectrum only is not reliable.

Given the large number of MWDs that have been identified since that study, and the mass availability of precise trigonometric parallax measurements for basically every object in the solar neighborhood \citep{Gaia2016, Gaia2018, Gaia2021}, it is now more than necessary to revisit/determine the atmospheric parameters of all known MWD in a homogeneous fashion as such an analysis is bound to provide important clues as to the nature and evolution of magnetic white dwarfs in general.

As a first step, this paper will look at all the objects that have been classified as magnetic hydrogen-rich white dwarfs in the Montreal White Dwarf Database \citep{MWDD}.
In \autoref{s:observations} we present the observations and our sample selection and in \autoref{s:theory} we describe our theoretical framework and the description of our analysis method for the photometric as well as the spectroscopic part of this work.
\autoref{s:results} present the details of our spectrophotometric analysis while \autoref{s:discussion} discusses those results and implications on our understanding of the evolution of MWDs.

\section{Observations} \label{s:observations}

\subsection{Sample Selection} \label{ss:sample}

As a first step, we aim to analyze in a homogeneous fashion all known magnetic hydrogen-rich white dwarfs that have been identified in the literature over the years \citep[][and references therein]{Ferrario2015a, Ferrario2020}.
We also include the many objects that have been flagged as DAH by visual inspection of thousands of SDSS spectra \citep[mostly from][]{Kepler2015, Kepler2016}.
Most of these magnetic white dwarfs have never, to our knowledge, been the subject of a detailed analysis.
Furthermore, we restrict ourselves to objects showing clear signs of magnetic line splitting ($B$ greater than about 2 MG), as the analysis of weakly MWDs requires either much higher resolution and signal-to-noise ratio spectroscopic observations than what is currently available (the bulk of our sample being from SDSS) or spectropolarimetric measurements \citep[such as][and references therein]{Vennes2018, Kawka2019, Bagnulo2021, Bagnulo2022}.
We thus first selected all the objects with the spectral type DA and \verb|Magnetic| flag in the Montreal White Dwarf Database \citep[hereafter MWDD,][194 stars]{MWDD} or with spectral type DAH (510 objects).
Stars that did not clearly show signs of magnetic splitting, or were known to have a magnetic field strength well below 1 MG based on published values in the literature were removed from our sample (see below).
We also removed a few objects that were modeled as a DA+DC unresolved binary \citep{Rolland2015} as the analysis of such systems requires special treatment.
While it is possible that a few genuine high field magnetic DA white dwarfs may have been missed, given that the MWDD includes all the MWD reported in the compilations of \citet{Ferrario2015a, Ferrario2020} as well as the new identifications from \citet{Kepler2015, Kepler2016}, we are confident that our sample contains almost all known hydrogen-rich MWD with fields higher than 2 MG that have been identified in the literature.

This leaves us with a sample of 661 objects.
However, a few of these objects had no optical spectra available to us, leaving us with a final sample of 651 potentially MWDs to examine.
Each of these objects is then analyzed in a homogeneous fashion using the best available spectra (see MWDD), astrometric measurements from Gaia DR2/EDR3 \citep{Gaia2016, Gaia2018, Gaia2021} and photometric measurements (either SDSS \emph{ugriz} \citep{York2000}, Pan-STARRS \emph{grizy} \citep{PanSTARRS} or \emph{JHK} and \emph{BVRI} \citep{Bergeron1997, Bergeron2001}, see \autoref{ss:photo-spectro-data}).

\autoref{f:color_mag_gaia} shows an observational Hertzsprung-Russell (H-R) diagram for the objects in our sample (colored dots) as well as WDs within 100 pc (black dots).
Evolutionary sequences for non-magnetic DAs are also shown.
The classification of the colored dots is the result of our analysis described in \autoref{s:results}.
The figure shows a clear separation between the MWD DAs and the rejected MWD candidates (probably normal DAs, see below), with MWDs tending to have lower luminosity, smaller radius, and therefore a higher mass than their non-magnetic counterparts.

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/reclassified_color_mag_gaia}
    \caption{\label{f:color_mag_gaia} Gaia H-R diagram for our sample of magnetic DAs. Black dots are the white dwarfs from the Gaia EDR3 within 100 pc. The various color dots represent our full sample of 651 stars, classified into categories explained in \autoref{s:results}.}
\end{figure}

\subsection{Photometric and Spectroscopic Data} \label{ss:photo-spectro-data}

In order to perform an analysis as homogeneous as possible, we decided to use SDSS \emph{ugriz} point-spread function magnitudes in priority when available.
When \emph{ugriz} data were not available, we relied on Pan-STARRS \emph{grizy} or previously published \emph{JHK} and \emph{BVRI} photometry, in that order of priority.
In the end, out of 651 objects, only 29 (4.5\%) had no \emph{ugriz} data.
Hence we used \emph{grizy} data for 26 stars and \emph{BVRI} photometry for the last 3.

The majority (598 out of 651, 92\%) of the spectra used for this work are collected from SDSS.
The spectra of the 53 objects without SDSS spectroscopy were taken from \citet{Limoges2013, Limoges2015}, \citet{Bergeron1997, Bergeron2001}, \citet{Gianninas2011}, \citet{Giammichele2012}, \citet{Tremblay2020} or archival data obtained by the Montreal group over the years.

\subsection{Note on Naming Convention} \label{ss:names}

Numerous names are often used to design a given white dwarf. In order to avoid confusion, we decided to follow a convention similar to that of \citet{Coutu2019} and use names based on ICRS coordinates at epoch and equinox 2000 (alternative names for any given star can easily be obtained through the search function on the Montreal White Dwarf Database).
Stars are thus referenced with JHHMM$\pm$DDMM, where HHMM corresponds to the right ascension (R.A.) in hours and minute and DDMM to the declination (decl.) in degrees and minutes in sexagesimal notation.
\autoref{t:names} provides the correspondence between the names of all objects in our sample with Gaia source ID, MWDD ID, and the coordinates in decimal format (degrees).

\section{Theoretical Framework} \label{s:theory}

\subsection{Magnetic Synthetic Spectra Calculations}

The presence of a strong magnetic field has a strong effect on both the atomic energy levels, line strengths as well as the radiative transfer as polarization has to be taken into account.

We rely on data provided by one of us \citep[S.Jordan, similar from][]{Forester1984, Wunner1985, Wunner1989} which provide transitional energies, wavelengths, and various quantum quantities needed to compute oscillator strengths for a wide range of magnetic field intensities from 0 to over 9000 MG.
The data used in this work were never fully published, and are denser than the ones available in the literature.
As can be observed in \autoref{f:balmer}, line positions vary substantially in a non-linear way for fields above about 10 MG.
The spectral features for objects with such strong fields can thus become much more complex than the usual Zeeman triplet found at lower fields.
We implemented an interpolating scheme developed by \citet{Steffen1990} to obtain the component data for arbitrary field intensities, with a linear extrapolation when required.
This extrapolation occurs at over 9000 MG for most components, although some weaker components only have calculations up to 1000 MG.
As most stars in our sample have a dipolar magnetic field $B_p$ of under 100 MG with only a very limited number of stars with $B_p$ over 500 MG, this extrapolation is not expected to be frequent.
We also obtained similar data for helium transitions and implemented them into our routines.
This work assumes only pure hydrogen atmospheres (other compositions will be explored in Hardy et al., in preparation).

An appropriate polarized radiative transfer routine also had to be implemented for a correct description of MWD radiation field.
\citet{Jordan2003} describe four methods to numerically solve the radiative transfer equations, and we decided to use their MATEXP method (Section 2.4 of their work), as it yielded smoother polarization spectra, was more numerically stable for our case, and was one of the fastest in computation time.
While in principle our code is now suitable for the analysis of other Stokes parameters, we restrict our analysis to the Stokes parameter $I$ (the luminous intensity) in this work as no polarimetric data is yet available for the majority of our sample\footnote{It should be noted that this is not the same as \emph{neglecting} the other Stokes parameters, as observational data is not yet widely available.}.

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/MAGLINH}
    \caption{\label{f:balmer} Position of the Balmer lines between 3000 and 8000 \AA{} as a function of magnetic field strength.}
\end{figure}

As a first-order approximation, we model the surface magnetic field of the stars with an inclined and offset dipole:

\begin{equation}\label{e:dipole}
    \mathbf{B}=\frac{3\left(\mathbf{m}\cdot \hat{\mathbf{r}}\right)\hat{\mathbf{r}}-\mathbf{m}}{{r}^{3}}
\end{equation}

\noindent where:

\begin{itemize}
    \item[] $\mathbf{B}$ is the magnetic field vector, in MG,
    \item[] $\mathbf{m} = \frac{B}{2} \hat{\mathbf{m}}$ is the dipole moment vector ($B$ is the dipole intensity, in MG, $\hat{\mathbf{m}}$ the unit vector of the dipole moment),
    \item[] $\hat{\mathbf{r}}$ is the unit vector of the position vector $\mathbf{r}$, from the dipole to where the field is measured (here, the surface of the star),
    \item[] $r$ is the absolute value of the position vector $\mathbf{r}$, in stellar radii.
\end{itemize}
We orient the magnetic moment according to the inclination, and apply the dipole offset by shifting the position vector accordingly.
The parameters used to describe the dipolar geometry (\autoref{f:dipole_schema}) are the dipole intensity $B_p$, which is the magnetic field intensity at the pole of a centered dipole, at the stellar surface ($1 \, R_*$ from the origin), the inclination angle $i$ between the dipole moment $\vec{m}$ and the observer's line of sight, and the dipole offset $a_z$ (in stellar radii $R_*$), the distance between the center of the dipole and the center of the star.

We discretize the stellar surface into many surface elements.
The number of surface element has to be large enough to correctly sample the variation of the magnetic field on the surface, but not too large as to slow down the calculation significantly.
We find that using 149 surface elements was optimal as no significant change to the computed synthetic spectra were found in tests using a finer resolution.
Once the specific intensity for each surface element has been calculated, we integrate across the stellar surface to obtain the Eddington flux for any given geometry.

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/dipole_angles_schema}
    \caption{\label{f:dipole_schema}Schema for the various parameters describing our dipole geometry. In red are the magnetic field vectors at the stellar surface, with the angle $\psi$ in green and $\theta$ in blue used to model each surface elements. $\vec{m}$, the inclination angle $i$ and dipole offset $a_z$ describe the geometry of the field.}
\end{figure}

In order to recover the magnetic field geometry for every star in our sample, a significant number of synthetic spectra across a wide range of effective temperatures, surface gravity, and geometries were required.
Pre-generating all the necessary dipolar-geometry intensity spectra would have been much more time consuming (and less expandable in the future) than integrating them on-the-fly using a precomputed grid of spectra with a constant field.
We therefore built a spectrum grid across 5 parameters: effective temperature (5000 K to 30,000 K in 15 steps), surface gravity ($\log{g}$ from 7.0 to 9.5 in 6 steps), field intensity ($10^{-2}$ MG to $10^4$ MG in 159 steps), the angle between the line of sight and the normal of the surface ($\psi$ from 0 to $\pi$ in 6 steps), and the limb darkening ($\mu = \cos{\theta}$ from 0 to 1 in 4 steps).
From this precomputed grid (85 GB of data), a synthetic spectrum with any desired field strength and geometry (not limited to offset dipole) can be reconstructed rapidly.
Only offset dipole geometries are considered in this exploratory work.

\subsection{Photometric Analysis Technique} \label{ss:photo}

The first step in our analysis is to estimate the effective temperature and surface gravity of our objects.
However, the presence of a strong magnetic field can affect significantly the determination of these atmospheric parameters.
As a start, we first obtained effective temperature and surface gravity with models with no magnetic field.
The atmospheric structures and synthetic colors used are those of \citet{Bergeron2019} \citep[same atmospheric models used by][]{GenestBeaulieu2019}.
Note that those models do not include the effect of magnetism on the stellar structure.
In fact, it is not clear at the moment how exactly the structure is affected by the presence of a strong magnetic field and as a consequence, all past studies, to our knowledge, used non-magnetic structure to calculate magnetic synthetic spectra \citep[see][for example]{Kulebi2009}.
We thus employed a similar strategy here but improved on the fitting technique by including explicitly the effect of magnetism when determining atmospheric parameters from photometric data (see below).

We obtain the effective temperature and surface gravity with zero-field models using an evolutionary algorithm \citep[based on][]{Storn1997} to fully explore the parameter space of the grid (this is a test to make sure we recover the non-magnetic parameters obtained with standard technique before proceeding with magnetic synthetic spectra).
We applied de-reddening corrections to the observed magnitudes to take into account the amount of interstellar material on the line of sight \citep{Schlafly2011}.
We then constructed new synthetic grids with various magnetic field strengths to account for the shifting of lines from one photometric band to the other as the field increases.
The atmosphere structures are the same as those used previously but the synthetic spectrum used to derive the synthetic photometric colors include the effect of the magnetic field on the line splitting.
With these photometric grids in hand, we proceeded to obtain a value for the atmospheric parameters for every star in our sample for field strength of 0, 10, 50, 100, 500, 700, and 900 MG.
For each star in our sample, we built an interpolating table for the atmospheric parameters (effective temperature, surface gravity), plus mass and radius, so we can get interpolated values for arbitrary magnetic field intensities that will be obtained in the next step.
While this is done assuming a constant surface magnetic field, this method allows to at least capture the main effect of the presence of strong magnetic fields, that is the shifting of lines in adjacent bands (the exact details of the geometry is a second-order effect in that matter).

\autoref{f:magnitudes} shows how the synthetic \emph{ugriz} colors are affected by the inclusion of magnetic splitting for various field strength and effective temperatures.
In the case of magnetic DA stars, the \emph{u} and \emph{g} bands are the most affected, since the Balmer lines (except H$\alpha$) are near them.
With weak or no magnetic field, there are barely any Balmer lines in the \emph{u} band, but a magnetic field with an intensity above a few MG pushes many components inside the band.
A similar phenomenon happens for the \emph{g} band, which is in the range of the Balmer lines starting from H$\beta$.
Here the many split components affect the \emph{g} band considerably. With increased magnetic field, the \emph{r} band is the least affected as H$\alpha$ components tend to stay inside the band and the \emph{i} band can only have H$\alpha$ components drift inside.
The \emph{z} band is a little different, since some of the Paschen lines are not split in our model, due to a lack of data.

Including the effect of magnetism in the synthetic spectra can thus have a non-negligible effect on the atmospheric parameter determination (effective temperature and surface gravity).
This is the first time, to our knowledge, that this effect is explicitly taken into account for the determination of the atmospheric parameters.
Uncertainties on the fitted values were obtained by fitting the photometry again, but with a least-squares minimization method that also returns the estimated uncertainties.
We used the values obtained with the fitting method described above as starting values for the algorithm, ensuring uncertainties around those values.

\begin{figure}
    \centering
    \includegraphics[width=0.8\linewidth]{compare_magnitudes_field_8.0}
    \caption{\label{f:magnitudes}Relative differences in magnitude for models with and without magnetic field.}
\end{figure}

\subsection{Spectroscopic Analysis Technique} \label{ss:spectro}

When fitting a given star, the process to obtain a synthetic geometry-aware spectrum is the following:
\begin{itemize}
    \item Provide the dipole geometry parameters: $B_p$, $i$, $a_z$;
    \item Generate an $N$-element discretized surface geometry and compute the average surface field intensity;
    \item Interpolate the effective temperature and surface gravity to that average field value;
          \begin{itemize}
              \item See previous \autoref{ss:photo}.
          \end{itemize}
    \item Use those atmospheric parameters with each surface element to obtain the $N$ spectrum elements by interpolating in the 5-D grid;
    \item Sum the $N$ spectra into the dipolar-geometry-aware spectrum at the desired atmospheric parameters.
\end{itemize}

We use this process in the fitting technique to obtain every spectrum needed for each iteration until convergence.
Due to the high dimensionality of the problem and the highly degenerate nature of the solutions, standard minimization methods using the gradient descent are ill suited for this situation, since many local minima are found with different dipolar geometry.
The usual steepest-descent algorithm would get stuck on the first minimum it would encounter based on an initial guess value.
Making this initial guess value would be different for each star and is not optimal.
This is why the optimization algorithm used relies on a stochastic nature and allows searching for a wide area inside the parameter space.
It is part of the \verb|scipy| Python package, which is an implementation of the algorithm by \citet{Storn1997}.

We tested our method by generating non-magnetic and magnetic synthetic stars with various magnetic field intensities and geometries and fed it to our fitting code.
Our algorithm was able to retrieve the geometry for both the magnetic and the non-magnetic stars with satisfactory precision.
For the special case of non-magnetic synthetic stars, we still used the algorithm to get a dipolar geometry, but the dipolar intensity from the solution is minimal and the other parameters (angle and offset) are less constrained and are of no physical value.

Typical uncertainties were estimated by changing the dipole parameters ($B_p$, $i$, $a_z$) until the resulting synthetic spectrum diverged significantly from the observations, for a typical magnetic star in our sample.
We observed that the relative uncertainties for the different parameters were very different, meaning the solutions were more sensitive to an increase of magnetic field intensity ($B_p$) than the other geometry parameters.
We found a relative uncertainty of 5\% for the $B_p$ parameter, 20\% for the inclination $i$, and 50\% for the dipole offset $a_z$.

In summary, the fitting technique for magnetic spectra begins with the creation of a random set of ($B_p$, $i$, $a_z$) parameters.
For each individual in this population, we go through the process described above to obtain a spectrum, and calculate the $\chi^2$ value against the observations.

\section{Spectrophotometric Analysis} \label{s:results}

Using the fitting technique described in \autoref{s:theory}, we fitted all 651 objects in our sample in a homogeneous way assuming single DA synthetic spectra with a dipolar geometry.
Each solution was then carefully inspected by eye to judge the quality of the fit. As a result of this laborious exercise, we could classify all our objects in broadly 5 categories (see \autoref{t:class_count}): i) stars that are probably not magnetic (400) ii) stars well fitted with a dipole geometry (140) iii) stars that are clearly hydrogen-rich magnetic white dwarfs but that are moderately well modeled with a dipolar geometry (45) iv) stars that have hydrogen features too shallow (28) and finally v) stars that we were simply unable to fit using the current theoretical framework (38). We now discuss each category in turn.

\input{class_count}

The first category contains more than half of our original sample (400 stars).
Most of these object had been classified as DAH or DAH: from visual classification of low signal-to-noise ratio spectra from various data release catalog of spectroscopically confirmed white dwarf stars from the Sloan Digital Sky Survey \citep{Kleinman2004,Kleinman2013,Kepler2015,Kepler2019}.
While visual inspection indeed seems to show possible shifted components near the core of some Balmer lines (thus the DAH classification), it quickly became clear that those classifications did not hold up against detailed comparison with magnetic synthetic spectra.
For example, it was possible in a few cases to force a magnetic field strength/geometry to fit some glitches near the core of one Balmer line only to find that such configuration lead to predicted components that are not observed in the other lines.
When considering all the lines simultaneously, our algorithm was not able to find a magnetic solution that was statistically better than a non-magnetic fit. In the end, although it is quite possible that some of these stars will be found to be genuinely magnetic white dwarfs when higher signal-to-noise or spectropolarimetric observations become available, we believe that classifying these objects as DAH based on these low quality SDSS spectra was perhaps a bit premature.
We report in \autoref{t:B_lt_1MB} the photometric atmospheric parameters (effective temperature and $\log{g}$) assuming $B = 0$ and reassign them with the spectral type DA until proof to the contrary.

\input{B_lt_1MG_table}

The next category comprises 140 objects for which an off-centered dipole model did a pretty decent job at reproducing the various magnetically shifted components.
\autoref{f:good_fit} shows an example of such stars.
While in some cases some minor discrepancies could be observed between the model flux and observations, the overall quality of the fit is deemed sufficient to label them as \emph{Good dipole fit}.
On the other hand, the quality of the dipole fit was clearly inferior for the third category (45 objects) that we call \emph{Bad dipole fit}.
\autoref{f:ish_fit} shows an example of such bad fits.
That so many objects cannot be well reproduced with a dipolar geometry should not be surprising.
For one, a shifted dipolar field geometry is probably only a first-order approximation for a much more complex field structure and we should not expect all objects to be easily modeled this way.
Second, it is well known that many MWDs show rapid photometric variations on timescales of minutes to hours \citep[see for example][and references therein]{Williams2022, Kilic2021}.
A good example is that of G183$-$35 \citep{Kilic2019}, an unusual white dwarf that shows an H$_\alpha$ line split into five components, instead of the usual three components seen in strongly magnetic white dwarfs (we also find five components for J1328+5908 and J1430+2811). Time-resolved spectroscopy of this object seem to support the idea that we are witnessing rotational modulations of a complex magnetic field structure.
As our sample is composed mostly of SDSS spectra that were integrated on a period of time that span an unknown fraction of the total period, it is no surprise that many objects are not well reproduced with a simple dipole. These objects, listed in \autoref{t:results_bad}, represent good candidates to study the complexity of magnetic field structure through time-resolved spectroscopy.

Nevertheless, although the overall quality of the fits for these stars are clearly inferior in terms of reproducing the exact shape of all the Balmer components, the resemblance between the best fit models and the observations suggests that the determined mean surface fields are probably in the right ball park.
As a consequence, we can be confident that the determined effective temperatures and surface gravities/masses for these objets are also closer to reality than those determined from non-magnetic models. All our fits are presented in Appendix \hyperref[f:all_fits]{A1}.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{"figures/fits/_VV2006__J120150.2+614256"}
    \caption{\label{f:good_fit} An example of a star that is successfully reproduced by a dipolar field geometry. The left panel is our best fit to the Balmer lines. Upper right panel represents the best fit to the $ugriz$ photometry while the bottom right panel shows the distribution of magnetic field elements at the visible surface of the star for our best fit solution.}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{figures/fits/SDSS_J084906.22+003722.8}
    \caption{\label{f:ish_fit} An example for a star where a shifted dipole provides a bad fit.}
\end{figure}

Another interesting category found within our sample contains objects with hydrogen features that are very shallow compared to what magnetic model spectra predict for the fitted atmospheric parameters and amount of splitting observed (28 objects, \autoref{t:shallow}).
\autoref{f:shallow_features} shows an example of such shallow features.
It must be noted that no combination of inclination, offset or field strength can reproduce those spectra.
It is also highly unlikely that a more complex multipole expansion could alleviate the issue as the depth of the central component of each line has only a weak dependence on the geometry.
In the case of J2257+0755 (EGGR156), \citet{Kulebi2009} managed to obtain a good fit by increasing the effective temperature to 40,000 K.
This solution, however, is clearly incompatible with the shape of the spectral energy distribution ($ugriz$) which requires an effective temperature in the vicinity of 14,000 K.
Lowering the effective temperature below $\sim 10,000$ K also reduces the depth of the absorption features, but again, at the cost of not reproducing the $ugriz$ photometry anymore.
Similar objects with shallow features have been modeled has DA+DC unresolved binary in \citet{Rolland2015} but in the case of J2257+0755, the observed photometric modulations rules out a close binary scenario and the star is rather interpreted as the rapidly rotating remnant of a double degenerate merger \citep{Williams2022}.
If so, rapid rotation may also produce an unusual surface composition that may explain the shallow lines found in such stars.
For example, an attractive solution to the absorption lines depth problem is that such objects have a mixed H-He composition.
We experimented with mixed atmosphere and indeed found that it was possible to reduce significantly the depth of the lines in that manner.
We calculated models with increasingly more helium, until the depth of hydrogen lines were well reproduced, and found that a 12,000K helium-dominated atmosphere with $\frac{N\left({\rm He}\right) }{N\left({\rm H}\right)} = 10^4$ can reproduce H$\alpha$ while not predicting helium features that are not observed.
Unfortunately, such models badly reproduce H$\beta$.
It is possible that the atmosphere is not homogeneous, with patches of different composition and temperature, but exploring the full range of possibilities is beyond the scope of this paper.
Again, time-resolved spectroscopy may eventually help us uncover the real nature of such objects.
We will thus postpone the analysis of these stars to future studies and exclude them from our discussion below.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{figures/fits/EGGR_156}
    \caption{\label{f:shallow_features} J2257+0755 (EGGR 156), an example of a star whose observed features are too shallow, according to our models. Due to the fitting algorithm trying to minimise the numerical difference between the synthetic spectrum and observations with shallow features, the magnetic structure here is not representative of reality.}
\end{figure}

\input{shallow_features_table}

Finally, the last category contains stars that could not be fit within our theoretical framework (38, \autoref{t:nonDA}).
These stars are all probably magnetic white dwarfs of some sort since non-magnetic objects at their fitted effective temperature and surface gravity, as obtained from fitting the $ugriz$ photometry, would show strong hydrogen or helium lines that would be easily recognized.
An example of such a star is presented in \autoref{f:magnetic_cant_fit} where none of the features could be matched using hydrogen-rich models with offset dipole.
The algorithm converged in the parameter space at a place where the surface elements field strength are spread over a large range, effectively washing out spectral features over the entire optical range.
While some of the multitude of features produced this way may match those in the observed spectrum, most of it is probably just random coincidences.
We thus cannot confirm the hydrogen-rich nature of these stars based on our minimum chi-square solution.
It is possible that our failure to find a good solution is due to some large temporal variations of the surface field strength with rotation over timescale much shorter than the integration time of the SDSS spectra.
Combined with that, it is also possible the geometry used here (inclined offset dipole) is too simple an assumption and a much more complex geometry is required.
Alternatively, the main atmospheric composition of those objects might simply not be hydrogen.
In the next paper of that series (Hardy et al. in preparation), we will explore helium-rich compositions and reassess the nature of these objects.

\input{magnetic_non_DA_table}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{figures/fits/2MASS_J02114816+2115491}
    \caption{\label{f:magnetic_cant_fit} An example of a magnetic white dwarf which we were unable to find a suitable dipolar geometry.}
\end{figure}

\input{results_good_table}
\input{results_bad_table}

\section{Discussion} \label{s:discussion}

From our original sample, we are left with 185 objects that we could fit with various degrees of quality using a dipolar geometry.
\autoref{f:compare_teff} shows that including the effect of the magnetic field in the synthetic color calculations tends to systematically increase the determined effective temperature compared to that obtained from non-magnetic models (the difference can be as high as 3700K for stars with the largest field intensity).
As explained in \autoref{ss:photo}, this is a direct consequence of the shifting of absorption lines across adjacent photometric bands as the magnetic field strength increases.
To compensate for this increase in effective temperature, it becomes necessary to decrease the stellar radii in order to reproduce the photometric data.
Hence, the derived masses are also increased by up to \Msol{0.1} compared to masses determined from non-magnetic models (see \autoref{f:compare_mass}).

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/compare_teff_magnetic}
    \caption{\label{f:compare_teff} Comparison of effective temperatures obtained from non-magnetic vs. magnetic synthetic spectra. The color scale indicate the strength of the magnetic field.}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/compare_mass_magnetic}
    \caption{\label{f:compare_mass} Comparison of stellar masses obtained from non-magnetic vs. magnetic synthetic spectra. The corresponding mass distribution histograms are also shown on the right and upper part of the figure. }
\end{figure}

It is interesting to discuss our sample in light of the results drawn from the 40 pc sample recently presented by \citet{Bagnulo2022}.
According to their interpretation, MWDs form two distinct populations distinguished by their very different typical masses.
The most massive magnetic white dwarfs, according to their hypothesis, would mostly be the result of binary mergers, producing strong magnetic fields that would be present at the surface as soon as the star enters the cooling sequence.
Such mergers are expected to be rapidly rotating.
On the other hand, young MWDs with mass below about $\sim$ \Msol{0.75} appear to be extremely rare and magnetism would gradually appear only after 2 Gyr or 3 Gyr, possibly due to the relaxation of a preexisting field buried in the interior or generated by a dynamo mechanism triggered by the onset of core crystallization \citep{Isern2017, Schreiber2021}.

\autoref{f:mass_dist} takes a closer look at the mass distribution of our sample as well as all the stars that we rejected in the previous section.
First, we see that the stars previously classified DAH that we dismissed as magnetic form a relatively broad distribution centered on $\sim$\Msol{0.6}, as expected for a population of non-magnetic white dwarfs.
This unusually large distribution can easily be explained by the fact that those stars are almost twice as distant, on average, than the stars that we successfully fit with a dipole.
These were thus objects with the very lowest signal-to-noise ratio spectra of the sample, and the observed splitting that led to their classification as DAH was most probably just noise.
Being more distant, the uncertainties on the photometric data and parallax measurements were also much larger, explaining the extremely broad mass distribution.

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/mass_distribution}
    \caption{\label{f:mass_dist} Mass distribution for the stars in our sample.}
\end{figure}

Next, we observe a sudden increase in the number of MWDs at a mass of about \Msol{0.7}.
This may indicate that the most massive MWDs have a different origin than the ones at lower mass, as suggested by \citet{Bagnulo2022}.
We must, however, remain cautious with this interpretation as our sample is mostly drawn from a magnitude-limited survey that certainly suffer from a number of selection biases difficult to quantify.
Nevertheless, a look at the distance distribution of our sample in \autoref{f:distance_mag_dist} suggests that the preponderance of massive white dwarfs in our magnetic sample is real and not merely an artifact of various selection effects.
Indeed, while it can be argued that massive white dwarfs may be, due to their much smaller radii, underrepresented in a magnitude limited sample, the fact that we see a cutoff in the number of magnetic stars with mass above \Msol{0.75} starting at about 250 pc suggest that we are perfectly capable of identifying such objects in SDSS at smaller distances.
Meanwhile, MWDs less massive than \Msol{0.75} are much larger (and thus brighter) and would be, at similar effective temperature, easily identify in that same volume, assuming that whatever selection bias affecting our sample is the same at all white dwarf masses.
According to the reasoning of \citet{Bagnulo2022} described above, this is totally to be expected if magnetism in lower mass objects only appears much later on the cooling sequence.
The big difference with our sample is that contrary to what is observed in the 40 pc sample, low mass MWDs are not completely absent for ages younger than 1-2 Gyr (see \autoref{f:mass_age}).
Taking only objects within 250 pc does not change the overall picture emerging from \autoref{f:mass_age} much; young MWDs less massive than \Msol{0.75}, while much rarer than heavier ones, do exist.
Their magnetic field strength, however, never reach values above $\sim$30 MG as is regularly found at higher masses.
It is possible, however, that whatever mechanism is responsible for the presence of strong magnetic fields in young low mass white dwarfs is sufficiently rare for it to not be statistically present in the 40 pc sample and that only looking at a much larger volume does this additional channel manifest its presence.
To summarize, our sample does suggest white dwarfs with strong magnetic field and mass above \Msol{0.75} are much more numerous than their lower-mass counterpart.
While selection biases may play a role in the exact proportion, this assessment is unlikely to change with a large volume-limited sample.
The young low mass MWDs in our sample could represent a separate channel not observed in the 40 pc sample that is superposed to the one where magnetism only start to show up after the onset of crystallization in the core.

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/distance_mag_distribution}
    \caption{\label{f:distance_mag_dist} Distance distribution of confirmed Magnetic DA stars, divided into two groups: One with masses below \Msol{0.75}, the other with masses above.}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/distribution_mass_age}
    \caption{\label{f:mass_age} Stellar mass vs. age, color coded with the mean surface magnetic field for the confirmed stars in our sample. The blue and orange lines represent the onset of crystallization for stars with a thick ($q_{He}=10^{-2}$, $q_H=10^{-4}$) and thin ($q_{He}=10^{-2}$, $q_H=10^{-10}$) hydrogen layer, respectively.}
\end{figure}

Are there other signatures for different channels/origin that emerge from the homogeneous analysis of the sample?
\autoref{f:mass_field} shows little correlation between the mean dipole field intensity (mean field intensity on the visible surface of the star) and stellar mass, except for the fact that all the largest fields are found only in high mass white dwarfs.
The stars that were badly fit with a dipole are slightly more prevalent among the most massive ($M \ge$ \Msol{0.75}) and most magnetic ($B \ge$ 10 MG) but can nevertheless be found at any mass, field or effective temperature.
We can see a slight tendency for heavier stars to be hotter, as there is a higher density of $> 25000$ K stars near \Msol{1.2}, although this could be a selection bias since more massive stars are smaller, and need to be hotter and more luminous to be seen at a given magnitude.
Looking instead at the effective temperature vs. the mean dipole field intensity (\autoref{f:teff_field}) does not highlight any discernable pattern that would allow us to distinguish between different groups that would possibly have different origins.
Such distinction may require a detailed analysis based on time-resolved spectroscopy in order to pinpoint possible differences in field structures.

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/distribution_mass_field}
    \caption{\label{f:mass_field} Mass vs. magnetic field intensity for the confirmed magnetic stars in our sample. The color scale shows effective temperature for each star.}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/distribution_teff_field}
    \caption{\label{f:teff_field} Effective temperature vs. magnetic field intensity for the confirmed magnetic stars in our sample.}
\end{figure}

\section{Conclusion} \label{s:conclu}

In this study, we used state-of-the-art magnetic synthetic models to perform the first large homogeneous analysis of hydrogen-rich MWDs since \citeyear{Kulebi2009}.
We carefully examined 651 white dwarfs that were labeled as magnetic or DAH compiled in the Montreal White Dwarf Database.
Since most of the sample is drawn from the Sloan Digital Sky Survey, only white dwarfs with magnetic fields large enough to produce an undeniable line splitting (above $\sim$1-2 MG) could be assessed.
From our original sample, 400 stars are reclassified as non-magnetic white dwarfs as the alleged splitting reported in the literature did not hold up against attentive scrutiny.
We find that 140 stars could be relatively well fitted with an offset dipole while another 45 were clearly DAH but a dipolar fit was passable.
We were also unable to fit 38 objects with our theoretical framework, suggesting that these objects either have a much more complex field structure or are simply not hydrogen-rich.
Finally, the depth of the absorption line for 28 stars were too shallow compared to model predictions, possibly due to non-homogeneous/mixed surface compositions.

We find a sharp increase in the number of MWDs starting at a mass of about \Msol{0.70}.
The strongest magnetic fields are also found only among the most massive white dwarfs (no field above 30 MG are found for $M \leq$ \Msol{0.70}).
Stars that are well and badly fitted with a dipole geometry appear to be disseminated evenly across the whole range of masses, effective temperatures and magnetic field strengths, leaving little clues about their possibly different origins.
The stars badly fitted with a dipole could have a much more complex field structure or simply vary on short time scaled compared to the spectra integration time.
Signature of the evolutionary origin of the magnetic field, if present, will have to wait for the availability of time-resolved spectroscopy in order to better quantify the details of their field structure.

A study of the 40 pc sample by \citet{Bagnulo2022} recently suggested that magnetism in white dwarf star with masses below \Msol{0.75} only start to appear after about 2 Gyr, possibly as the result of a dynamo mechanism that is triggered by the start of core crystallization.
However, based on our sample, it appears that there is an additional channel that is capable of producing very young magnetic white dwarfs with $M \leq$ \Msol{0.70}.
With about two dozen such object found in a volume much larger than the 40 pc sample, it may, however, not be statistically significant to not find any in the local sample.
Because our sample certainly suffer from numerous selection biases, we refrain from quantifying the proportion of such objects as a function of age or mass.

To conclude, detailed analysis of large samples of magnetic white dwarfs are needed to provide new insights on the nature and evolution of such stars.
Our study, while inconclusive, represent the largest effort in that endeavor since the work of \citet{Kulebi2009}, almost doubling the number of magnetic white dwarfs meticulously analyzed in a homogeneous fashion.
Our investigation established the existence of various families of magnetic objects that deserve further scrutiny, especially those not fit well with a dipole geometry or those with shallow absorption features.
Some object that may not have a hydrogen-rich surface composition will be examined in detail in the next paper of this series.
Time-resolved spectroscopy, spectropolarimetry and much higher signal-to-noise/resolution observations are highly awaited for advancing our understanding of the origin and evolution of magnetic white dwarfs.
The unique models and fitting tool developed for this project have certainly proven themselves as reliable, and will be the basis of studies that will push even further our understanding of MWDs.

\section*{Data Availability}

The data underlying this article are available on the Montreal White Dwarf Database\footnote{\url{https://www.montrealwhitedwarfdatabase.org/}} \citep{MWDD}.

\bibliographystyle{apj}