\section{Introduction} \label{intro}

Type Ia supernovae (SNe Ia) are the thermonuclear explosions of carbon-oxygen (C/O) white dwarfs (WDs) \citep{Iben1984, Nomoto1982}. SNe Ia have been used to calibrate cosmological distance scales \citep{Nomoto1997, Perlmutter1999} and cosmological models \citep{Riess1998}, yet our knowledge on the progenitors of SNe Ia remains incomplete. Although the precise progenitor system is unknown, it is generally accepted a primary WD is stimulated via mass accretion from, or merges with, a companion to explode as a SN Ia \citep{Wang2012reviewofprogenitors}.

In the so-called ``single-degenerate'' (SD) scenario \citep{WhelanIben1973, Nomoto1982}, the companion is a non-degenerate hydrogen- or helium-star (H-star and He-star, respectively) which donates material from its outer layers to the primary, mass-accreting WD. In the ``double-degenerate'' (DD) scenario \citep{Iben1984, Webbink1984}, the donating companion is another WD. A merger is possible in either scenario. Whether or not a primary WD and companion star will lead to the SD or DD scenario depends on parameters such as the initial masses of the components, their separation, and metallicity \citep{Iben1994, Wang2010metallicity}. These parameters also determine which of the many sub-channels each scenario will likely take \citep{Wang2011, Liu2018}.

The explosion of the primary WD can be triggered multiple ways in either scenario \citep{Iben1984, Webbink1984, Pakmor2012, manganese, Papish2015, Shen2018a}. Theory and observations both suggest different combinations of triggering mechanisms may describe various observed abundance patterns in SN Ia remnants \citep{manganese, Fink2014}. This implies that both near-Chandrasekhar (near-Ch) mass WDs and sub-Chandrasekhar (sub-Ch) mass WDs, either in the SD or DD scenario, explode to contribute to all observed SNe Ia.

Recent models suggest that the donor will almost always survive the explosion \citep{Pan2012, Shen2018b, Liu2018}. In cases of complete primary WD detonation, the surviving companion will be kicked away with a velocity roughly equal to the binary's pre-explosion orbital velocity \citep{Eldridge2011, Wang2011}. This population of kicked non-degenerate/degenerate stars should show evidence of the SN Ia in their photosphere \citep{Pan2012, Pan2012b, ShenSchwab2017, Shen2018b, Tanikawa2018}. Regardless of the evolutionary stage of the kicked donor, it will ultimately evolve to the main WD cooling sequence \citep{Hansen2003, Justham2009}. Therefore a population of kicked WDs are expected from the SN Ia scenario, harbouring evidence of their close proximity to thermonuclear supernova. We denote this entire kicked WD population as donor remnants (DRs) regardless of their evolutionary phase during the progenitor stage.

It has also been suggested that the primary WD can survive its own failed detonation in a sub-luminous SN Iax event \citep{Foley2013, Kromer2013}. If the primary WD remains intact, the asymmetry in the partial deflagration could lead to its ejection from the binary and tell-tale ashes may remain visible in its photosphere \citep{Jordan2012}. It has also been argued that the partial deflagration does not warrant an ejection of the companion \citep{Kromer2013, Fink2014}. In either case, the possibly intact primary WD is a different kind of remnant -- one that should show evidence of the sub-luminous SN Iax in their own photosphere. These ``primary remnants'' (PRs) should also eventually fall to the main WD cooling sequence \citep{Jordan2012, Kromer2013, Zhangetal2019}.

LP~93-21\footnote{Also known as LHS~291 \citep{Luyten1976}, EGGR~434 \citep{Greenstein1977}, WD~1042+593 and WD~1042+59 \citep{McCook1999}, and SDSS~J104559.13+590448.3, SDSS~J104559.14+590448.3, and SDSS~J104559.15+590448.2 \citep{SDSSdr14}.} is a single, DQ spectral-type WD, known since 1976 to have strong carbon Swan Bands and a peculiarly high proper motion \citep{Luyten1976}. In this paper, we examine the properties of LP~93-21 and argue that it is an evolved PR which was kicked to its abnormally high space velocity by its own partial deflagration. We also discuss the possibility that LP~93-21 may also be an evolved donor. In Section \ref{Observations}, we present recent observational data and atmospheric analysis of LP~93-21 to provide updated estimates of its total space motion, mass, age, and origin. In Section \ref{discussion}, we analyse the kinematics of LP~93-21 in context with these updated estimates and assess the viability of the SNe Ia/Iax progenitor scenario. We conclude in Section \ref{Conclusion}.

\section{Observations} \label{Observations}

Here we report on the most up-to-date observations and parameter estimates of LP~93-21. LP~93-21 was first examined in detail by \citet{Greenstein1977} as a carbon degenerate. Its spectrum showed particularly strong C$_{2}$ molecular bands, and its proper motion made it a high-velocity WD. The Sloan Digital Sky Survey (SDSS) has spectroscopic confirmation of over nineteen-thousand white dwarfs, most with reported radial velocities \citep{Kleinman2013, SDSSoverview, SDSSdr14}. The European Space Agency's \textit{Gaia} satellite has also observed over twenty-five thousand white dwarfs, measuring their proper motion and parallax with milliarcsecond accuracy \citep{Gaiamission, Lindegren2018, DR2}. Both surveys have since re-observed LP~93-21, providing improved measurements which confirm both its peculiarly fast motion and strong carbon spectral features.

\subsection{Kinematics} \label{kinematics}

The high precision of \textit{Gaia} DR2 measurements (see Table \ref{measuredproperties}) for LP~93-21 provide a parallax measurement of $\varpi = 17.48 \pm 0.14$ mas that can be directly converted to distance $d = 57.2 _{-0.3} ^{+0.4}$ pc without significant loss of accuracy \citep{BailerJones2018, Luri2018}. We may calculate the 3-D space motion for LP~93-21 by combining \textit{Gaia}'s proper motion measurements and the SDSS radial velocity measurement. We include the radial velocity measurement from its SDSS spectrum of $462 \text{ km}\text{ s}^{-1}$ \citep{Kleinman2013} and add a conservative radial velocity error of $\pm$ $20 \text{ km}\text{ s}^{-1}$, which we assume is uncorrelated with the \textit{Gaia} observables.

To integrate LP~93-21's orbit, we used the \texttt{astropy} \citep{astropy2013, astropy} affiliated \texttt{python} package \texttt{gala}\footnote{\url{https://gala-astro.readthedocs.io/en/latest/index.html}} \citep{Price-Whelan2017}. We chose the Milky Way potential in \texttt{gala} which assumes a mass-model containing a spherical nucleus and bulge based on the Hernquist potential for a spheroid \citep{bulge}, a disk based on the Miyamoto-Nagai potential for a flattened mass distribution \citep{Disk}, and a spherical NFW dark matter halo potential \citep{NFW}. The relative weighting for each component is taken from (\cite{Bovy2015}; see their Table 1). We define the Sun's position at $x = -8.3$ kpc from the Galactic centre in the midplane ($z = 0$) of the Milky Way, and use an upper-limit circular velocity estimate of 250 km s$^{-1}$ \citep{Bovyetal2012, Eilersetal2019}. We drew initial positions for each orbit from the observed covariance matrix, as described in \citet{Lindegren2018}, which factors in the correlation between the \textit{Gaia } observables.

We find LP~93-21 is currently travelling at approximately $v_{\textrm{gal}} = 605 \text{ km}\text{ s}^{-1}$ with a local galactic escape speed of approximately $v_{\textrm{esc}} = 557 \text{ km}\text{ s}^{-1}$, making it locally unbound to the Milky Way's gravitational potential. At a distance of $d = 57.2 _{-0.3} ^{+0.4}$ pc, LP~93-21 is the closest hyper-runaway star to the Sun. Figure \ref{orbit} shows LP~93-21's integrated orbit forwards and backwards in time 250\,Myr -- approximately the timescale of one Solar orbit around the Galactic centre. For all 100 orbits drawn, LP~93-21 is unbound to the Milky Way's gravitational potential. Figure \ref{orbit_zoom} shows its proximity to the Sun. LP~93-21's closest approach to the Sun would have been approximately 70,000 years ago at a distance of about 46 pc.

The observed velocity of LP~93-21 is currently being influenced by the Milky Way's gravitational potential. By integrating LP~93-21's orbit in \texttt{gala}, we take into account its velocity over time only as it has been affected by the Milky Way's gravitational potential.
If we consider an extragalactic origin for LP~93-21 then as it approached the disk of the Milky Way the gravitational pull began to accelerate it and increase its speed (see Figure \ref{evolution}, bottom right).

If LP~93-21 did have an extragalactic origin, then by integrating backwards we find that
at a distance of approximately 100 kpc from the Galactic centre (taken as a proxy for the virial radius of the Milky Way), LP~93-21's galactocentric velocity is about 400\,km s$^{-1}$. This would correspond to a flight time of about 220\,Myr from a 100 kpc galactocentric radius to LP~93-21's current position. This provides an upper limit for ejection time if the SN kick occurred within the Milky Way.

\begin{table*}
\begin{tabular}{@{}lcccc@{}}
\toprule
\multicolumn{1}{c}{Parameter} & Symbol & Value & Units & Source \\ \midrule
Right Ascension & $\alpha$ & $161.49 \pm 0.08$ & deg & \citet{DR2} \\
Declination & $\delta$ & $59.07 \pm 0.10$ & deg & \citet{DR2} \\
Proper motion in Right Ascension & $\mu_{\alpha}$ & $-1019.19 \pm 0.14$ & mas yr$^{-1}$ & \citet{DR2} \\
Proper motion in Declination & $\mu_{\delta}$ & $-1462.53 \pm 0.17$ & mas yr$^{-1}$ & \citet{DR2} \\
Parallax & $\varpi$ & $17.48 \pm 0.14$ & mas & \citet{DR2} \\
Radial velocity & $v_{\textrm{rad}}$ & $462 \pm 20$ & km s$^{-1}$ & \citet{SDSSdr14} \\
Apparent $G$-band magnitude & $G_{\textrm{app}}$ & 17.7 & mag & \citet{DR2} \\
Absolute $G$-band magnitude & $G_{\textrm{abs}}$ & 13.8 & mag & \citet{DR2} \\
$BP-RP$ colour & $C_{\textrm{bp-rp}}$ & 0.24 & mag & \citet{DR2} \\
Mass & $M$ & $1.029 \pm 0.015$ & M$_{\odot}$ & \citet{Kilic2018} \\
Effective temperature & $T_{\textrm{eff}}$ & $8690 \pm 120$ & K & \citet{Kilic2018} \\
Surface gravity & $\log \textrm{g}$ & $8.701 \pm 0.02$ & cm s$^{-2}$ & \citet{Kilic2018} \\
Carbon abundance & [C/He] & $-3.51$ & dex & \citet{Kilic2018} \\
White dwarf cooling age & $t_{\textrm{WD}}$ & $2.715 \pm 0.08$ & Gyr & \citet{Kilic2018} \\ \bottomrule
\end{tabular}
\caption{Measured properties of LP~93-21.}
\label{measuredproperties}
\end{table*}

\begin{figure*}
    \centering
    \includegraphics[width=1.0\textwidth]{LP_93_21_855361055035055104.png}
    \caption{LP~93-21's orbit. The current position of LP~93-21 is shown by a black circle. The integrated orbits trace LP~93-21's movement forward in time (red) and backwards in time (blue). The spread in orbit directions is caused by the uncertainty in \textit{Gaia } and SDSS measurements. The extent of the Milky Way's spiral arms are approximated by the darker, dashed circle, and the galactic centre is approximated by the intersection of the lighter, dashed lines. The direction to the Large Magellanic Cloud is denoted by a black arrow. Each hash mark indicates approximately 25 million years of flight either forwards or backwards in time. Note LP~93-21 is unbound in all 100 integrated orbits.}
    \label{orbit}
\end{figure*}

\begin{figure*}
    \centering
    \includegraphics[width=1.0\textwidth]{LP_93_21_zoom_855361055035055104.png}
    \caption{LP~93-21's proximity to the Sun. The current position of LP~93-21 is shown by a black circle. The position of the Sun is denoted by a yellow star. The integrated orbits trace LP~93-21's movement forward in time (red) and backwards in time (blue). The spread in orbit directions is caused by the uncertainty in \textit{Gaia } and SDSS measurements. LP~93-21's closest approach around 70,000 years ago was approximately 46 pc.}
    \label{orbit_zoom}
\end{figure*}

Based on kinematic information alone, we can exclude the Galactic centre as the origin of LP~93-21, as well as the Magellanic clouds, and six known supernova remnants (SN1006A, SN1054A, SN1604A, SN185A, SN393A, and Vela; \cite{Green2014}).
Figure \ref{evolution} shows the kinematic and positional evolution for LP~93-21. The galactocentric z-coordinate $Z_{\textrm{gal}}$, galactocentric radius $R_{\textrm{gal}}$, and heliocentric distance $ d_{\textrm{helio}} $ all show asymptotic behavior, characteristic of unbound stars.

\begin{figure*}
    \centering
    \includegraphics[width=0.7\textwidth]{LP_93-21_evolution.png}
    \caption{Kinematic and positional parameters for LP~93-21 shown as a function of flight time $t$. Movement forward in time is denoted in red and movement backwards in time is denoted in blue. Top left: galactocentric z-coordinate $Z_{\textrm{gal}}$. The vertical extent of the thick disk is approximated by the dashed lines. Top right: galactocentric radius $R_{\textrm{gal}}$. Bottom left: heliocentric distance $d_{\textrm{helio}}$. Bottom right: galactocentric rest-frame velocity $v_{\textrm{gal}}$.}
    \label{evolution}
\end{figure*}

\subsection{Current mass \& age estimates} \label{currentestimates}

\citet{Kilic2018} identified 142 \textit{Gaia } sources as halo WDs and presented detailed model atmosphere analysis for each one based on \textit{Gaia } parallaxes and optical and near-infrared photometry. They used a pure helium atmosphere model with trace amounts of carbon to simultaneously fit the spectral energy distribution (SED) and the optical spectra for LP~93-21. The results of this fit are summarised in Table \ref{measuredproperties}. \citet{Weidemann2005} notes the inaccuracies of using He-dominated DB WD models as a replacement for DQ WDs without accounting for this enhanced carbon composition, implying that the detailed analysis by \citet{Kilic2018} is the current state-of-the-art. LP~93-21's mass is found to be 1.03\,$\text{M}_{\odot}$, making it a massive outlier in SDSS (\cite{Kepler2007}; see their Figure 13). LP~93-21's WD age is estimated by \citet{Kilic2018} to be 2.715 Gyr, under the assumption of zero mass transfer.

\citet{Leggettetal2018} presented new US Naval Observatory Flagstaff Station parallaxes for over 170 WDs. They combined recent \textit{Gaia } parallaxes with photometry spanning from the mid-infrared to the ultraviolet to determine flux-calibrated SEDs for each WD. The SEDs were compared to flux distributions provided by various recent model atmospheres compiled by the authors. \citet{Leggettetal2018} demonstrate the models reproduce the full SED very well for the entire sample of WDs. LP~93-21's mass is again found to be an outlier at 1.14\,$\text{M}_{\odot}$, and its WD age is estimated to be 2.366 Gyr, under the assumption of zero mass transfer.

\begin{figure*}
    \centering
    \includegraphics[width=0.7\textwidth]{HRD_updated.png}
    \caption{Hertzsprung-Russel diagram showing a random sample of one million \textit{Gaia } stars with usual features such as the main sequence, giant branch, and white dwarf sequence. LP~93-21 is denoted by a black dot within the main white dwarf sequence at an approximate absolute \textit{G}-band magnitude of 13.8 mag. Other known unbound WDs are shown. The three ``D$^6$'' hyper-runaway WD candidates \citep{Shen2018b} are denoted by blue dots. The hyper-runaway WD candidate LP 40-365 \citep{Raddi2018} is denoted by a red dot. The hyper-runaway US 708 \citep{Hirsch2005} is denoted by a green dot.}
    \label{HRD}
\end{figure*}

\section{Discussion} \label{discussion}

LP~93-21 is an exceptionally fast-moving WD, moving much quicker than would be expected for a normal, Milky Way halo WD \citep{Oppenheimer2001, Bergeron2003haloWDs, Bergeron2005highvelocityhaloWDs, Pauli2006WDkinematics, Ducourant2007haloWDs}. LP~93-21's observed velocity likely excludes it from normal dynamical ejection from a globular cluster because simulations estimate that $\approx$ 99\% of all dynamical ejections of this type have less than $200~$km s$^{-1}$ velocities \citep{Perets2012}. LP~93-21's observed velocity also likely excludes it from the core-collapse binary supernova ejection scenario. Simulations estimate that only the fastest few percent of stars ejected via core-collapse supernovae reach speeds of over $200$ km s$^{-1}$ \citep{tauris2015}.

If LP~93-21 originates from the Milky Way, we find its velocity was very likely never below 400 km s$^{-1}$. Therefore, the mechanisms that may explain LP~93-21's velocity (if from the Milky Way) include dynamical ejection from the Galactic centre or a binary SNe Ia/Iax ejection. Given LP~93-21 does not originate from the Galactic centre (or any other location where a central, massive black hole is expected to be, see Figure \ref{orbit}), and can be excluded from other dynamical scenarios based on its observed speed, the most plausible mechanism to describe its velocity is the SNe Ia/Iax ejection scenario.

The progenitor configuration most compatible with LP~93-21's observed speed is the SD scenario proposed for SNe Ia/Iax between a primary WD and He-burning donor star \citep{Wang2009, Wang2017heliumAccretion, wongetal2019, Baueretal2019}. In the SN Ia scenario, the primary WD reaches Chandrasekhar-mass through accretion, explodes as a SN Ia, and the DR is ejected at roughly the pre-SN orbital velocity \citep{Eldridge2011, tauris2015, Shen2018b}. In the SN Iax scenario, ejection speeds on the order of a few tens to a few hundred km s$^{-1}$ are expected for both the DR and PR, due to the expected asymmetric explosion and ejecta interaction \citep{Jordan2012, Pan2012impact, Pan2012evolution, Liu2012ms, Liu2013heliumCompanions, Kromer2015}.

It is unlikely that LP~93-21 was a WD while donating mass to its former primary companion in the DD scenario. Due to the relatively short upper-limit Milky Way flight time of 220 Myr for LP~93-21 (see Section \ref{kinematics}), its maximum Milky Way speed was likely never above $\approx\,625$km s$^{-1}$ (see Figure \ref{evolution}). The expected tell-tale kinematic signature (i.e. $>$ 1,000 km s$^{-1}$ in all cases) indicates LP~93-21 should still be possessing a much quicker total speed if it were ejected as a WD in the DD scenario \citep{Shen2018b, Tanikawa2018}. This is also indicated by Figure \ref{orbit}. If shot through the Milky Way's disk with a characteristic speed of $>$ 1,000 km s$^{-1}$, LP~93-21 should have only been accelerated to quicker speeds after its ejection. Thus, the DD scenario is unlikely for LP~93-21 based on its observed speed alone.

The SD SN Ia scenario implies LP~93-21 was a mass-donating He-star companion that was ejected by the detonation of its primary\footnote{In the same way, the SN Iax scenario also implies LP~93-21 may have been a former He-star donor.}. In this scenario the relatively short upper-limit Milky Way flight time of 220 Myr for LP~93-21 (see Section \ref{kinematics}), when combined with the updated WD cooling age estimates of over 2\,Gyr from both \citet{Kilic2018} and \citet{Leggettetal2018}, demands an origin from elsewhere than the Milky Way. LP~93-21 would have then evolved into a WD while travelling through inter-galactic space.

We find this scenario is an unlikely explanation for LP~93-21's origin because of the low probability that an extra-galactic hyper-runaway star would pass so close to the Sun.
We took over 850 nearby galaxies from the Updated Nearby Galaxy Catalog \citep{Karachentsev2013galaxycatalog} and estimated SN Ia rates given the galaxy mass and rate-size relations from \citet{Lietal2011}.
We assume an average speed of 600\,km\,s$^{-1}$, which assumes some deceleration as the DR escapes the gravitational potential of the host galaxy and acceleration as the DR approaches the Milky Way. If we assume that the DRs are ejected in isotropic directions and that a fully-cooled down DR could be observed by \textit{Gaia } out to a distance of 100~pc, then the expected number of  observable DRs from a the $n$th galaxy is given by the ratio of observable volume to the volume of a sphere with a radius equal to the galaxy's distance
\begin{equation}
    E_n = \left(\frac{100\,\textrm{pc}}{d_n}\right)^3 \textrm{SNuM}(M_n)\Delta{t}
\end{equation}
\noindent{}where $d_n$ is the distance to the host galaxy, and $\textrm{SNuM}(M_n)$ is the SN Ia rate as a function of galaxy mass. Here $\Delta{t}$ is conservatively taken as the maximum time the DR could take to pass through the observable volume: the time taken for DR to travel twice the observable radius (towards and away from us; $\Delta{t} = 2\frac{100\,\textrm{pc}}{600\,\textrm{km\,s}^{-1}} \approx 3.3\times10^6\,\textrm{yr}$).
After summing the expected number of observable DRs from all nearby galaxies we find that only $10^{-4}$ of fully-cooled, extra-galactic DRs would be expected. This value remains much less than one even if we relax our assumptions.
This low expectation value indicates it is highly unlikely any extra-galactic DR is within 100 pc of the Sun, including LP~93-21.

If ejected as a He-star, a corresponding rotation speed and detailed elemental abundances would be paramount in determining the viability of this scenario for LP~93-21. It is expected that LP~93-21 would be a very fast rotator, and that a nontrivial abundance of decayed nickel might still be visible on its surface \citep{Fuller2012, Pan2012impact, Liu2013rotation}. Trace amounts of intermediate-mass elements might also point to the He-star SD SN Ia channel as the explanation for LP~93-21's observed speed \citep{Pan2012evolution, Liu2013heliumCompanions}.

Conversely, the SN Iax explanation implies LP~93-21 was possibly a former \textit{primary} WD in the same SD configuration between a WD and companion He-star. The upper-limit Milky Way flight time is now consistent with its WD cooling age, by this explanation. Primary WDs which have possibly survived their own partial deflagrations make a new class of stars based on distinct kinematic and abundance signatures \citep{Raddietal2019}. The observed speeds of these objects agree well with LP~93-21's. But as shown by the red dot in Figure \ref{HRD}, these objects do not appear as ordinary WDs. The handful of objects found appear to be ``puffed-up'' WDs, meaning they are brighter and slightly red-shifted when compared to normal WDs (see Figure \ref{HRD} and \cite{Raddietal2019}). The degeneracy in their cores may have been lifted slightly due to the deflagration and accretion of ejecta \citep{Shen2018b, Raddietal2019}.

Simulations of WDs which include radiative levitation and gravitational settling suggest that the observable effects from the deflagration should dissipate over relatively short timescales for the primary WD, but that a second heating and brightening phase could last on the order of 10\,Myr until these objects eventually evolve back to the canonical WD cooling sequence \citep{Zhangetal2019}. The current location of LP 40-365 is therefore likely temporary, as the direct visual effects dissipate much sooner than a typical Galactic crossing time for an unbound PR (about 100$-$150 Myr, \cite{Raddi2018}).

Due to its position in Figure \ref{HRD}, LP~93-21 could represent the first fully ``cooled-down'' PR observed in this new class of stars similar to LP~40-365 \citep{Zhangetal2019, Baueretal2019}. LP~93-21's enhanced C abundance agrees well with this scenario, although its bulk composition may not fully reflect in surface abundances \citep{Zhangetal2019}. Carbon burning is required in Iax deflagrations, but incomplete in such explosions, so PRs are likely to have substantial carbon abundances \citep{Zhangetal2019}.  LP~93-21 is a massive WD outlier ($>$\,1.0\,$M_{\odot}$), and three-dimensional hydrodynamic simulations suggest that pure deflagrations of hybrid C/O/Ne WDs may leave massive bound remnants, which is inconsistent with the observed stars like LP~40-365 ($\approx$\,0.2\,-\,0.3\,$M_{\odot}$) which were likely left by pure deflagrations of C/O or O/Ne WDs \citep{Kromer2015, Raddietal2019}.

\section{Conclusions} \label{Conclusion}

Using \textit{Gaia } DR2, we were able to confirm LP~93-21 exhibits a total space velocity of approximately $v_{\text{gal}} = 605 \text{ km}\text{ s}^{-1}$, with a local galactic escape speed of approximately $v_{\text{esc}} = 557 \text{ km}\text{ s}^{-1}$. We found its orbit to be unbound to the Milky Way in all cases. LP~93-21's orbit implies that it has travelled at most 220 Myr if originating from the Milky Way. The SN Iax-ejection scenario is consistent with its observed velocity, and a detailed spectrum may confirm if LP~93-21's own partial deflagration in a SN Iax caused such an ejection from the Milky Way.

If a portion of SNe Ia/Iax progenitors are a result of the He-star donor channel, the post-explosion kinematics necessarily produces a population of single WDs with high velocities, possibly containing distinguishing elemental abundances \citep{Wang2009, Liu2018, Wang2018review, Zhangetal2019, Raddietal2019}. Thus it is possible that many high-velocity, slightly evolved He-stars and high-velocity WD remnants originate from both the SN Ia and SN Iax progenitor channels.

However, WDs are intrinsically very faint, so magnitude-limited surveys have historically had difficulties finding representative populations for WD analysis outside the Solar neighbourhood. LP~93-21's chance proximity to the Sun enabled discovery and observations easier than most for the likely numerous WD remnants kicked in the SNe Ia/Iax scenarios. In the future, magnitude-limited surveys should improve for very faint magnitudes and many more ``cooled-down'' remnants of these scenarios are likely to be found.

Regardless of the origin, LP~93-21 is the closest hyper-runaway WD candidate to the Sun. Other SNe Ia/Iax progenitor remnant candidates are shown on a Hertzsprung-Russel diagram of BP-RP colour versus absolute \textit{G}-band magnitude in Figure \ref{HRD}. In four cases, a detailed spectral analysis was used in conjunction with a kinematic analysis to constrain the nature of the progenitor systems that likely caused their high velocities. All four candidates were WDs found to be enhanced in key elements and unbound to the Milky Way's gravitational potential. Due to the scarcity of known objects of this type, it is possible that LP~93-21 is the first identified example of a fast-moving, ``cooled down'' primary remnant and could represent a missing piece of the SNe Ia/Iax-ejected WD cooling sequence.

Further spectral analysis would provide constraints on the viability of the SNe Ia/Iax scenario in explaining LP~93-21's observed speed. A detailed spectrum would allow measurements of LP~93-21's spin and of key elemental abundances that might or might not still be on its surface. More surviving remnants are likely contained within the \textit{Gaia } data and still need to be found throughout the Milky Way, waiting to help shed light on the SNe Ia/Iax progenitor problem.

\section*{Acknowledgements}

The authors would like to thank
    the anonymous referee for a detailed review,
    as well as
    Morgan Fraser (University College Dublin),
    Alexander Heger (Monash University),
    Adrian Price-Whelan (Princeton),
    Amanda Karakas (Monash University),
    Tin Long Sunny Wong (University of California Santa Cruz), and
    Paul Canton (University of Oklahoma)
for useful conversations. A.~R.~C. is supported by the Australian Research Council (ARC) through Discovery Project DP160100637.

This research has made use of NASA’s Astrophysics Data System (\url{https://ui.adsabs.harvard.edu/}). This research has made use of the following software: \texttt{astropy}, a community-developed core Python package for astronomy \citep{astropy2013, astropy}, \texttt{TOPCAT} \citep{topcat}, \texttt{gala} \citep{Price-Whelan2017}, \texttt{numpy} \citep{numpy}, and \texttt{matplotlib} \citep{matplotlib}.

This work has made use of CosmoHub \citep{cosmohub}. CosmoHub has been developed by the Port d’Informac\'io Cient\'ifica (PIC), maintained through a collaboration of the Institut de F\'isica d’Altes Energies (IFAE) and the Centro de Investigaciones Energ\'ieticas, Medioambientales y Tecnol\'ogicas (CIEMAT), and was partially funded by the “Plan Estatal de Investigac\'ion Cient\'ifica y T\'ecnica y de Innovac\'ion” program of the Spanish government.

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

This research has made use of data from the Sloan Digital Sky Survey (\url{https://www.sdss.org/}). Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges
support and resources from the Center for High-Performance Computing at
the University of Utah. The SDSS web site is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the
Participating Institutions of the SDSS Collaboration including the
Brazilian Participation Group, the Carnegie Institution for Science,
Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics,
Instituto de Astrof\'isica de Canarias, The Johns Hopkins University,
Kavli Institute for the Physics and Mathematics of the Universe (IPMU) /
University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory,
Leibniz Institut f\"ur Astrophysik Potsdam (AIP),
Max-Planck-Institut f\"ur Astronomie (MPIA Heidelberg),
Max-Planck-Institut f\"ur Astrophysik (MPA Garching),
Max-Planck-Institut f\"ur Extraterrestrische Physik (MPE),
National Astronomical Observatories of China, New Mexico State University,
New York University, University of Notre Dame,
Observat\'ario Nacional / MCTI, The Ohio State University,
Pennsylvania State University, Shanghai Astronomical Observatory,
United Kingdom Participation Group,
Universidad Nacional Aut\'onoma de M\'exico, University of Arizona,
University of Colorado Boulder, University of Oxford, University of Portsmouth,
University of Utah, University of Virginia, University of Washington, University of Wisconsin,
Vanderbilt University, and Yale University.

\bibliographystyle{mnras}

\section{Motivation for expected observables} \label{motivation}

In this section we present a short summary of the most likely observational signatures associated with a DR from a handful of the most likely SNe Ia progenitor scenarios. We also present an estimate to the number of DRs we would expect in our observable volume within the Milky Way.

\subsection{Donor remnants to SNe Ia}

Tidbit here about how both SD and DD are thought to contribute to the rate of all SNe Ia and even SNe Iax??? **

\subsubsection{Double-denegerate remnants}

Although the exact progenitors to SNe Ia remain uncertain, a sub-Chandrasekhar (sub-Ch) mass C/O WD accreting helium-rich mass from a WD companion in the DD scenario is one possible pathway (see \cite{Maoz2014} for a review). There is now evidence to suggest all DD WD systems eventually undergo dynamically \textit{unstable} mass transfer, even those that eventually coalesce to form a single merger remnant \citep{Pakmor2013, Shen2015, Danetal2015}. Unstable mass transfer may trigger a helium shell detonation in the primary WD and a subsequent detonation in its core during the initial phases of mass transfer, leading to a SN Ia well before complete tidal disruption of the companion WD \citep{ShenMoore2014, ShenSchwab2017}. This is the so-called ``double-detonation’’ channel thought to occur in the large majority of DD SNe Ia (also known as the ``D$^6$'' scenario or ``helium-ignited violent merger model’’, see \cite{Pakmor2013, Shen2018b, Shen2018a}). Dynamically \textit{stable} mass transfer to the primary WD is also possible \citep{Bildstenetal2007, ShenBildsten2009}, and can lead to a helium shell detonation with subsequent core detonation in a variant of the double-detonation scenario.

The double-detonation channel implies the companion WD survives the explosion of the primary in either stability regime, leaving a surviving DR traveling at high velocity \citep{Shen2015, ShenSchwab2017, Tanikawa2018}. This is supported by observational constraints. If the companion WD were to fully tidally disrupt, it should form a torus surrounding the primary WD \citep{Danetal2014}. This torus should cause significant asymmetry to the SN Ia ejecta -- which is directly at odds with the high levels of symmetry found in SNe Ia remnants \citep{Bulla2016}. The material in the torus surrounding the primary WD should also pollute the SN Ia ejecta. This material should be observable by emission lines but is almost never detected \citep{Shen2018b}.

The double-detonation channel also implies different effects (over different timescales) for the DR in either SNe Ia mass transfer stability regime. For the unstable regime, shock heating from the SN Ia and its ejecta, tidal field heating, and the capture, accretion, and delayed radioactive decay of $^{56}$Ni $\longrightarrow$ $^{56}$Co $\longrightarrow$ $^{56}$Fe likely affect the outer layers of the companion WD \citep{ShenSchwab2017, Tanikawa2018}. These effects manifest by increasing the luminosity of the WD DR (i.e. $>$\,0.1 L$_{\odot}$) and causing the colour to become bluer than expected anywhere from a few years to a few centuries \citep{Iben1998, ShenSchwab2017, Shen2018b}. By extension, these effects are therefore important mostly in discriminating younger DRs from ``local’’, historical SNe Ia rather than older DRs with ambiguous origins. For the stable regime, the donor is only mildly distorted by these effects and should return to the WD cooling sequence much sooner.

\citet{Tanikawa2018} simulate the interaction between, and the evolution of, a primary 1.0 M$_{\odot}$ C/O WD with a large He shell and a secondary 0.6 M$_{\odot}$ C/O WD. The He shell on the primary WD was artificially thick so that the He and subsequent C/O detonation would succeed easily. \citet{Tanikawa2018} followed the evolution of the binary system from the initiation of the He detonation to the effects of the SN Ia ejecta on the expected DR observables. The DR's mass remains roughly equal, but it does capture a small amount of material originating from the centre of the detonation in the primary WD (the slowest ejected material, most easily caught by the companion, should originate from the densest region). This material consists mostly of $^{56}$Ni and He because when a C/O detonation propagates through a high-density region it mostly synthesises $^{56}$Ni, and a small amount of He is left as a residual due to photo-dissociation. \citet{Tanikawa2018} found the DR WD contaminated with approximately 0.03 M$_{\odot}$ of $^{56}$Ni and 3\,$\times$\,10$^{-4}$ M$_{\odot}$ He, which is similar to the results of \citet{ShenSchwab2017}, the only other cited simulation following the D$^6$ SN Ia scenario and the subsequent ejecta impact on DR observables.

The $^{56}$Ni will undergo delayed radioactive decay $^{56}$Ni $\longrightarrow$ $^{56}$Co $\longrightarrow$ $^{56}$Fe within the DR \citep{ShenSchwab2017}. The decaying $^{56}$Ni can be identified as anomalous, but it is still unclear if the material stays on the surface of the DR due to radiative levitation \citep{Chayer1995a, Chayer1995b} or falls to the centre due to sedimentation \citep{Paquette1986, Dupuis1992}. It is necessary to perform further simulations to follow the time evolution of the DR to know whether or not decay products remain on the surface of the WD, and if so,en how long they might impact late-time observations.

\citet{Shen2018b} model the runaway velocities of DRs produced by the D$^6$ scenario. \citet{Shen2018b} use the results of \citet{ShenSchwab2017} and \citet{Shen2018a} to assume the runaway velocity of the donating WD (both He and C/O WDs in their simulations) at the point of SN Ia is roughly equal to the orbital velocity of the companion star during Roche-lobe overflow. It is found the runaway velocities of surviving D$^6$ DRs are remarkable, in every case more than 1,000 km s$^{-1}$. For the case of a 0.6\,$M_{\odot}$ companion WD to a 1.0\,$M_{\odot}$ primary WD, the DR runaway velocity is found to be 1800 km s$^{-1}$. For the case of a near-equal mass binary of a 1.0\,$M_{\odot}$ companion and a 1.0\,$M_{\odot}$ primary, the DR runaway velocity is found to be an astounding 2,200 km s$^{-1}$.

The visual changes from the SN and its ejecta should last on the order of tens to hundred of years for local DD DRs originating from the Milky Way. It is still unknown if the captured SN ejecta is easily observable on the surface of the DR at late-time observations. Therefore the most distinguishing property of a DR originating from either a Milky Way or extragalactic SN Ia in the DD scenario is its exceptional space velocity.

\subsubsection{Single-denegerate remnants}

The SD scenario is another possible SNe Ia progenitor channel, in which a sub-Ch mass C/O WD accretes either hydrogen- or helium-rich mass through Roche-lobe overflow, or through a strong wind, from the companion to reach Ch-mass and detonate as a SN Ia (see \citet{Maoz2014} for a review). The donating companion may be a main-sequence (MS) hydrogen burning star (H-star), naked helium star (He-star, often called a ``(hot) subdwarf'' due to its hydrogen envelope being stripped off by binary interaction), or red giant (RG). Different configurations for either secondary can lead to a detonation of the primary WD, but each case depends on the previous common envelope and mass transfer phases. The survival of an intact companion star in the SD scenario is also expected from theoretical simulations and observational constraints \citep{}. In the same way as before, the intact companion from the SD scenario should display evidence of the explosion in its atmosphere \citep{} and be a kinematic outlier \citep{}.

Regardless of the sub-channel, the main SD channel predicts the primary sub-Ch mass WD will accrete matter from the companion star and increase in mass to approach the Chandreskhar-mass limit ($\approx$1.44\,$M_{\odot}$). Once it reaches this mass limit, the primary reaches sufficient density at the core to initiate a thermonuclear runaway explosion. This neatly explains the uniformity of supernovae remnant (SNR) light curves observed \citep{}. But for MS H-star companions, only a fairly narrow range in accretion rates above 10$^{-7}\,M_{\odot}$ yr$^{-1}$ are allowed for stable hydrogen mass transfer and burning on the surface of the primary WD, which allow the primary to avoid a regular nova explosion \citep{}. This narrow window heavily restricts possible H-star companion configurations leading to SN Ia.

Observed DRs possibly originating from the H-star companion SD channel are still in question, whereas a few promising DD channel DRs have been presented, providing further evidence for the difficulty H-stars have in producing SNe Ia \citep{}. SNe Ia SNRs are also characterised by a lack of hydrogen in their spectra, which is naturally explained by the DD channel. But in the hydrogen-accretion model of the SD scenario, simulations suggest the impact of the SN may strip off $\sim 0.03-0.15$ $M_{\odot}$ of H-rich material from H-star companions and almost entirely removes the envelop of RG companions, polluting the SN ejecta with hydrogen \citep{}. This is directly at odds with observations as well. High-resolution spectroscopy of two relatively hydrogen-rich extraglactic SNe Ia estimated an upper limit of only 0.01\,$M_{\odot}$ possible for the amount of stripped hydrogen polluting the ejecta \citep{}.

In addition to the difficulty of H-star companions leading to SN Ia, and observational constraints indicating SNe Ia are unlikely byproducts of H-rich mass transfer, simulations also suggest that the SN Ia effects on the companion H-star would also make it difficult to detect large visual differences at late-time observations (as evolved DRs).
**Tidbit here about how the ejecta will affect the luminosity and visuals of the MS H-star on short-ish timescales
In the long pre-explosion mass transfer phase associated with H-rich accretion, the companion is pushed out of thermal equilibrium. The increased radius and luminosity of the companion arising from mass transfer in a binary system is known as ``thermal bloating'' \citep{Knigge2011}, but it is unlikely this effect lasts longer than the timescales associated with captured ejecta decay.
Still unknown  about the radiative levitation/settling - need more simulations!

**Although spin may be used as another diagnostic, apart from abundances. Jury's still out.

**The kinematics associated with H-star DRs also likely makes it difficult to detect them apart from background stars. --> Not fast enough to make them out on tangential velocity, which should be ``easy'' in big tangential velocity surveys. Not really fast enough in radial velocity alone either, and RV measurements are much less prevalent.

Following the difficulty of hydrogen-mass accretion in the SD scenario and the likely small, relatively difficult to distinguish (by abundances or kinematics) population of DRs associated with the hydrogen-mass accretion model, more attention in the past decade has been placed on the \textit{helium}-accretion model of the SD scenario \citep{}. It has been found that the He-star companion model to the SD channel can possibly describes the discrepancy in theory and observations for both SD and DD SNe Ia \citep{}. It has also been found to be a likely dominant pathway to the formation of (massive) DD systems \citep{}.

**He-mass transfer is easier for WDs \citep{}, hydrogen envelope can easily be stripped off via binary interaction \citep{}, the naked He-star gets closer in separation for easier He-mass transfer, different configurations of WD+He-star were studied and found XX amount lead to SN Ia themselves, also lots evolved to WD+WD and found XX amount lead to SN Ia as DD scenario \citep{}. When including the WD+He-star configuration to the DD scenario in this one particular advanced simulation, found to make up 2/3 of observed SNe Ia rate \citep{}.

The WD and He-star ``intermediate'' channel also neatly explains the historical delay-time-distribution (DTD) issues stemming from a singular SD or singular DD prescription. Delay time is defined as the time interval from star formation to the thermonuclear runaway explosion (thus a distribution for a population). DTDs predicted by theoretical works not including the WD and He-star channel have large deficits associated with the early (e.g. $\lesssim$ 1 Gyr) and old (e.g. $\gtrsim$ 8 Gyr) epochs of delay time \citep{}. Upon including the WD and He-star pathway, simulations show the DTDs of SNe Ia match better with observations, even more so at the old and early epochs \citep{}.

**Explain the WD+He-star scenario leading to SN Ia as an example, and explain why the helium mass being accreted is more efficient for SN Ia.

**Kinematics of this separate He-star kicked population versus the H-star population. Why distinctly different? Observations support this? (not really) \citet{Tauris2015, Yisikandeer2016} Possible reasons for discrepancy? What would we be looking for, then??

**Unfortunately, simulations also suggest the SN Ia ejecta effects on the ejected companion He-star would also make it difficult to detect large visual differences at late-time observations (as DRs). Examples of accreted ejecta material amounts here with how it would effect observations over short/long time scales.

\subsubsection{Sub-luminous SNe Ia remnants}

A peculiar sub-class of SNe Ia likely produced via the SD scenario have recently been presented in the literature, characterised (among other things) by fainter light curves, lower SN expansion velocities, and lower ejecta mass when compared to Branch-normal Ia events \citep{}. The lower energetics and lower mass ejecta associated with these events, dubbed SNe Iax, are thought to be explained by the partial or full deflagration of a primary WD in the SD scenario (depending on the energy release, \cite{}). This ``failed detonation'' (FD) scenario implies a partial deflagration can leave the primary WD intact and the conditions after the deflagration do not lead to a subsequent detonation (for energies e.g. < ). See \citet{} and \citet{} for discussion on avoiding the transition to detonation. At sufficient enough energies (e.g. > ), the FD scenario implies a full deflagration will propagrate through the entire primary WD and unbind it completely.

What nucleosynthesis is associated with partial/full deflagrations in this scenario? How will that effect the primary? Companion? What might we expect to see here

What kinematics are associated with partial/full deflagrations in this scenario? How might it affect the primary? Companion? What might we expect to see here

Examples in the literature of DRs from this scenario?

\subsection{Observable donor remnant population}


\section{Introduction}
\label{intro}

Mass transport in porous media is one of the central topics in hydrology. For example, mass transport has been intensively studied in the contexts of carbon dioxide capture and storage (\citealt{Benson_2008}) and enhanced oil recovery (\citealt{Thomas_2008}). Moreover migration of the soil pollution of radioactive cesium in Fukushima, Japan, which is considered as mass transport in porous media, is an urgent issue (\citealt{Matsuda_2015}). The phenomenon of mass transport in porous media has been analyzed by the advection-diffusion equation (ADE) since 1950s. However, different examples in which ADE is not applicable have been found. Adams and Gelhar showed that the mass transport predicted by ADE was significantly different from measured values in a field scale experiment (\citealt{Adams_1992}). In the laboratory scale, Berkowitz, Scher, and Silliman found that the transport in a laboratory flow cell does not obey ADE, but it shows anomalous diffusion (or non-Fickian transport) and deviates from Fick's law (\citealt{Berkowitz_2000}). In such flow of anomalous diffusion, unlike the prediction by ADE, the peak of concentration barely moves and its distribution has a long tail. The deviation from ADE was investigated by a random particle simulation (\citealt{Kennedy_2001}).

To reproduce anomalous diffusion, continuous time random walk (CTRW) has been developed since the first proposal by Montroll and Weiss (\citealt{Montroll_1965}). Anomalous diffusion in porous media in different situations can be reproduced by CTRW when waiting time and jump distance are suitably adjusted. Although CTRW is widely used in the study of porous media (\citealt{Berkowitz_1998,Hatano_1998,Levy_2003,Nissan_2017}), the relation between the movement of tracer particles and the waiting time and jump distance of CTRW is not apparent. By taking the correlation in waiting time into account, correlated continuous time random walk was proposed (\citealt{Montero_2007}). In addition to CTRW, the following models were proposed. In tempered anomalous diffusion, strong anomalous diffusivity is suppressed by the introduction of an exponential cutoff for waiting times (\citealt{Meerschaert_2008}). The stochastic hydrology was developed (\citealt{Gelhar_1986,Rubin_2003}). The mobile/immobile (MIM) model and fractional MIM model divide mass transport into the moving part and immovable part (\citealt{Genuchten_1976,Schumer_2003}). It is known that diffusion equations with fractional derivatives are obtained in the asymptotic limit of CTRW (\citealt{Metzler_2000}). In fractional advection-dispersion equation,  time and spatial derivatives are replaced by fractional derivatives (\citealt{Benson_2000}). The use of fractional derivatives has attracted attention as a theoretical tool for anomalous diffusion. Fractional derivative models are compared (\citealt{Wei_2016,Sun_2017}), parameters in fractional equations were estimated (\citealt{Chakraborty_2009,Kelly_2017}), the physical meaning of fractional derivatives was discussed (\citealt{Liang_2019}), and orders of fractional derivatives were taken to be variables (\citealt{Sun_2009}).

To study transport in porous media, in addition to the concentration in mass transport, actual flow was visualized (\citealt{Moroni_2009,Sen_2012}) and the flow was reproduced by large scale simulation (\citealt{Anna_2013,Liu_2016}). These early studies show that tracer particles and fluid particles change velocity as a function of time when along the flow path they make a detour around obstacles in the the porous medium.

Although CTRW and fractional equations have achieved certain success, the physical meaning of the random walk process and physical background of fractional-order derivatives are not easy to understand. In this paper, we propose to use the linear Boltzmann equation (LBE), in which the physical scattering process is more apparent. Figure \ref{fig:schematic} illustrates the Boltzmann transport picture of the transport in porous media. Tracer particles in a porous medium flow along fissures and fractures. We regard the process of a tracer particle changing its direction as a scattering. Thus for the ensemble of tracer particles, the transport in the porous medium can be viewed as the transport of noninteracting particles which are scattered by scatterers in different directions with a certain probability distribution.

\begin{figure}
\centering
\includegraphics[clip,width=0.4\textwidth]{Fig.1.eps}
\caption{
(Color online) The transport of tracer particles in a porous medium. This transport is regarded as noninteracting particles which undergo scatterings by scatterers.
}
\label{fig:schematic}
\end{figure}

The rest of the paper is organized as follows. In Sec.~\ref{sec:Method}, the linear Boltzmann equation is introduced. We explore the advection-diffusion approximation in Sec.~\ref{sec:ADA}. In Sec.~\ref{sec:ADO}, we solve the linear Boltzmann equation. Section \ref{sec:column} is devoted to column experiments. Conclusions are given in Sec.~\ref{sec:concl}.

\section{Linear Boltzmann transport}
\label{sec:Method}

The linear Boltzmann equation, radiative transport equation, or transport equation governs the transport of noninteracting particles with scattering and absorption. The equation describes neutron transport in a reactor (\citealt{Case_1967}) and light transport in random media such as clouds and the intersteller medium (\citealt{Apresyan_1996,Chandrasekhar_1960,Ishimaru_1978,Sobolev_1976,Thomas_1999}). Let us regard the transport of tracer particles along flow paths as the transport of particles in a homogeneous random medium. See the schematic figure in Fig.~\ref{fig:schematic}. There is a random arrangement of flow paths in beads and sand in column experiments. At a branching point where several paths split into different directions, tracer particles are {\em scattered} to these directions with certain probabilities. The average distance that a tracer particle travels along a flow path between branching points is an effective mean free path $v_0/\sigma_s$ (see (\ref{proposed_model}) for $v_0,\sigma_s$). Indeed, the use of the linear Boltzmann equation for the flow of tracer particles in porous media was proposed by Williams (\citealt{Williams_1992a,Williams_1992b,Williams_1993a,Williams_1993b}).

A remark is necessary for our Boltzmann transport. The (nonlinear) Boltzmann equation can be simulated by the lattice Boltzmann scheme (\citealt{Succi01}), and the lattice Boltzmann method has become one of commonly used numerical methods to investigate the flow in porous media. Instead of fluid in porous media, we focus on the Boltzmann transport of tracer particles in the fluid. Hence, the proposed equation (\ref{proposed_model}) below is a linear equation.

Having in mind the transport in column experiments which is described in Sec.~\ref{sec:column}, we consider the one-dimensional linear Boltzmann equation. Let $\mu\in[-1,1]$ be the cosine of the polar angle and $v_0>0$ be the inherent particle speed. Although the coefficient of the spatial derivative is $v_0\mu$ in Williams' study (\citealt{Williams_1992a,Williams_1992b}) , in this paper, we take into account advection, which will be denoted by $u\ge0$. Hence the velocity $v$ is written as $v=u+v_0\mu$. Let $\psi(x,\mu,t)$ be the angularly dependent number density of particles at position $x\ge0$ in direction $\mu$ at time $t\ge0$. We can write the linear Boltzmann equation as follows.
\begin{equation}
\label{proposed_model}
\frac{\partial}{\partial t}\psi(x,\mu,t)+(u+v_0\mu)\frac{\partial}{\partial x}\psi+(\sigma_a+\sigma_s)\psi=
\frac{\sigma_s}{2}\int^{1}_{-1}\psi(x,\mu,t)\,\mathrm{d}\mu.
\end{equation}
where isotropic scattering is assumed in the integral on the right-hand side. Here, $\sigma_a\ge0$ and $\sigma_s>0$ are the absorption and scattering coefficients, respectively. Let $n_0$ be the initial particle number density. Tracer particles are injected to the column at $x=0$ in the normal direction. The initial and boundary conditions are given by
\begin{eqnarray}
\label{analytical solution1_1}
&&\psi(x,\mu,0)=0,\\
\label{analytical solution1_2}
&&\psi(0,\mu,t)=n_0\delta(\mu-1),\quad\mu>-\eta,
\end{eqnarray}
where $\delta(\mu-1)$ is the Dirac delta function and
\begin{equation}
\label{analytical solution1_8}
\eta=\frac{u}{v_0}.
\end{equation}
We note that unlike the usual half-range boundary condition of $\mu>0$ (\citealt{Case_1967}), the boundary value in (\ref{analytical solution1_2}) is specified for $\mu\in(-\eta,1]$ due to the presence of the advection. We regard the column as the half space ($0<x<\infty$) and have $\psi\to0$ as $x\to\infty$.

The particle number density $n(x,t)$ is calculated as
\begin{equation}
\label{defnumdens}
n(x,t)=\int_{-1}^1\psi(x,\mu,t)\,\mathrm{d}\mu.
\end{equation}
Then the concentration is given by
\begin{equation}
C(x,t)=\alpha n(x,t),
\end{equation}
where $\alpha$ is the mass per particle.

\section{Advection diffusion approximation}
\label{sec:ADA}

Let us derive the advection-diffusion equation from the linear Boltzmann equation. According the usual procedure of the diffusion approximation (\citealt{Duderstadt_1979,Ishimaru_1978}), we assume that the angular dependence of $\psi(x,\mu,t)$ is weak and written as
\begin{equation}
\label{LBE_ADE_1}
\psi(x,\mu,t)=\frac{1}{2\alpha}C(x,t)+\frac{3}{2\alpha v_0}J(x,t)\mu,
\end{equation}
where,
\begin{align}
\label{LBE_ADE_2}
C(x,t)&=
\alpha\int^{1}_{-1}\psi(x,\mu,t)\,\mathrm{d}\mu,
\\
\label{LBE_ADE_3}
J(x,t)&=
\alpha v_0\int^{1}_{-1}\mu\psi(x,\mu,t)\,\mathrm{d}\mu.
\end{align}
In (\ref{LBE_ADE_1}), we implicitly assumed a large propagation distance and long time. If absorption is strong, $\psi$ becomes zero before the form (\ref{LBE_ADE_1}) is achieved. Hence we assume $\sigma_a$ is small. Indeed in Sec.~\ref{sec:column}, we will see that $\sigma_a$ is negligibly small. Indeed, the approximation in (\ref{LBE_ADE_1}) is called the $P_1$ approximation (\citealt{Case_1967}). Below, we will proceed further and obtain the diffusion equation. We substitute the form (\ref{LBE_ADE_1}) for $\psi$ in (\ref{proposed_model}). By integrating the resulting equation over $\mu$, we obtain
\begin{equation}
\label{LBE_ADE_4}
\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+\frac{\partial J}{\partial x}+\sigma_aC=0.
\end{equation}
By integrating the equation over $\mu$ after multiplying $\mu$, we obtain
\begin{equation}
\label{LBE_ADE_5}
\frac{1}{v_0}\frac{\partial J}{\partial t}+\frac{u}{v_0}\frac{\partial J}{\partial x}+\frac{v_0}{3}\frac{\partial C}{\partial x}+
\frac{\sigma_a+\sigma_s}{v_0}J=0.
\end{equation}
We assume that the first term and second term of the above equation are small and can be dropped. Then we have
\begin{equation}
\label{LBE_ADE_6}
J=-D\frac{\partial C}{\partial x},
\end{equation}
where
\begin{align}
\label{LBE_ADE_8}
D&=\frac{v_0\ell^*}{3},
\\
\label{LBE_ADE_9}
\ell^*&=\frac{v_0}{\sigma_a+\sigma_s}\simeq\frac{v_0}{\sigma_s}.
\end{align}
Hence we obtain the following advection diffusion equation by substituting (\ref{LBE_ADE_6}) into (\ref{LBE_ADE_4}).
\begin{equation}
\label{LBE_ADE_7}
\frac{\partial C}{\partial t}-D\frac{\partial^2 C}{\partial x^2}+u\frac{\partial C}{\partial x}+\sigma_aC=0,
\end{equation}
For the photon diffusion equation, $\ell^*$ is known to be independent of the absorption coefficient (\citealt{Furutsu_1994}). The diffusion approximation holds when $x$ is sufficiently larger than the transport mean free path $\ell^*$.

The diffusion equation (\ref{LBE_ADE_7}) reduces to the usual advection-diffusion equation when $\sigma_a=0$. Since we found $\sigma_a\approx0$ in our column experiments, the solution by Ogata and Banks (\citealt{Ogata_1961}) is used in this paper when the concentration calculated from the linear Boltzmann equation is compared to the advection-diffusion equation. In Appendix \ref{wabsorption}, the solution to (\ref{LBE_ADE_7}) is presented.

Finally, we consider the boundary condition. Let us substitute (\ref{LBE_ADE_1}) for $\psi(x,\mu,t)$ in the boundary condition (\ref{analytical solution1_2}). By integrating both sides of the resulting equation over $\mu$ from $-\eta$ to $1$, we obtain
\begin{equation}
\label{bc2}
C(0,t)-\frac{3D\eta(1-\eta)}{2}\frac{D}{u}\frac{\partial C}{\partial x}(0,t)=\frac{2\alpha n_0}{1+\eta}.
\end{equation}
Thus we obtain the Robin boundary condition. In this paper, we neglect the flux term and write
\begin{equation}
\label{bc1}
C(0,t)=C_0,
\end{equation}
where $C_0>0$ denotes the boundary source. Since the diffusion approximation does not hold on the boundary, the actual proportionality constant between $\alpha n_0$ and $C_0$ is a fitting parameter. We introduce
\begin{equation}
\beta=\frac{\alpha n_0}{C_0}.
\end{equation}

\section{Concentration from the linear Boltzmann equation}
\label{sec:ADO}

Let us consider the Laplace transform.
\begin{equation}
\label{analytical_solution_1_5}
\hat{\psi}(x,\mu,p)=\int^{\infty}_0\mathrm{e}^{-pt}\psi(x,\mu,t)\,\mathrm{d}t.
\end{equation}
We introduce
\begin{equation}
\label{analytical solution1_6}
\mu_t=\frac{\sigma_a+\sigma_s+p}{v_0},\quad
\mu_s=\frac{\sigma_s}{v_0}.
\end{equation}
Then our transport equation is written as
\begin{equation}
\label{analytical_solution_1_4}
(\eta+\mu)\frac{\partial}{\partial x}\hat{\psi}(x,\mu,p)+\mu_t\hat{\psi}(x,\mu,p)=\frac{\mu_s}{2}\int^{1}_{-1}\hat{\psi}(x,\mu,p)\,\mathrm{d}\mu,
\end{equation}
with the boundary condition,
\begin{equation}
\label{analytical solution1_9}
\hat{\psi}(0,\mu,p)=\frac{n_0}{p}\delta(\mu-1),\quad\mu>-\eta.
\end{equation}

We will compute $\hat{\psi}$ with the analytical discrete ordinates method \mm{(ADO)} (\citealt{Siewert_1999,Barichello_2000,Barichello_2001}).
Unlike the standard ADO, the coefficient $\mu_t$ is a complex number since the variable $p$ is complex and moreover, the boundary condition is given for $\mu>-\eta$ instead of $\mu>0$.
We begin by decomposing $\hat{\psi}$ into the ballistic term $\hat{\psi}_b$ and scattering term $\hat{\psi}_s$ as
\begin{equation}
\label{analytical solution2_3}
\hat{\psi}(x,\mu,p)=\hat{\psi_b}(x,\mu,p)+\hat{\psi_s}(x,\mu,p).
\end{equation}
Here, $\hat{\psi}_b$ obeys
\begin{equation}
\label{analytical solution2_4}
(\eta+\mu)\frac{\partial}{\partial x}\hat{\psi_b}+\mu_t\hat{\psi_b}=0,
\quad x>0,\;-1\le\mu\le 1,
\end{equation}
with the boundary condition
\begin{equation}
\hat{\psi_b}(0,\mu,p)=\frac{n_0}{p}\delta(\mu-1),
\quad\mu>-\eta,
\end{equation}
and $\hat{\psi}_s$ obeys
\begin{equation}
\label{analytical solution2_7}
(\eta+\mu)\frac{\partial}{\partial x}\hat{\psi_s}+\mu_t\hat{\psi_s}=\frac{\mu_s}{2}\int^{1}_{-1}\hat{\psi_s}\,\mathrm{d}\mu+q,
\quad x>0,\;-1\le\mu\le1,
\end{equation}
with the boundary condition
\begin{equation}
\hat{\psi_s}(0,\mu,p)=0,
\quad\mu>-\eta,
\end{equation}
where
\begin{equation}
\label{analytical solution2_9}
q(x,\mu,p)=\frac{\mu_s}{2}\int^{1}_{-1}\hat{\psi_b}(x,\mu,p)\,\mathrm{d}\mu=
\frac{n_0\mu_s}{2p}\mathrm{e}^{-x\mu_t/(\eta+\mu)}.
\end{equation}
We obtain
\begin{equation}
\label{analytical solution2_6}
\hat{\psi_b}(x,\mu,p)=
\frac{n_0}{p}\mathrm{e}^{-x\mu_t/(\eta+\mu)}\delta(\mu-1).
\end{equation}

Let us drop $p$ and write $\hat{\psi}_s(x,\mu)=\hat{\psi}_s(x,\mu,p)$ and $q(x,\mu)=q(x,\mu,p)$. For the computation of $\hat{\psi}_s$, we discretize the integral by the Gauss-Legendre quadrature and obtain
\begin{align}
\label{analytical solution3_5}
&
(\eta+\mu_i)\frac{\partial}{\partial x}\hat{\psi_s}(x,\mu_i)+
\mu_t\hat{\psi_s}(x,\mu_i)
\nonumber \\
&=
\frac{\mu_s}{2}\sum^{N}_{j=1}w_j
\left[\hat{\psi_s}(x,\mu_j)+\hat{\psi_s}(x,-\mu_j)\right]+q(x,\mu_i),
\end{align}
where $\mu_i,w_i$ ($i=1,2,\dots,2N$) are abscissas and weights, respectively. These $\mu_i,w_i$ with $0<\mu_1<\cdots<\mu_N<1$ and $\mu_{N+i}=-\mu_i$ ($i=1,\dots,N$) can be calculated by the Golub-Welsch algorithm (\citealt{Golub_1969}). In the numerical calculation in Sec.~\ref{sec:column}, we set
\begin{equation}
N=30.
\end{equation}
Furthermore, we introduce $N_{\eta}$ as the largest integer such that $-\eta<\mu_{N_{\eta}}$. The scattering part $\hat{\psi}_s$ is obtained as
\begin{equation}
\label{analytical solution3_7}
\hat{\psi_s}(x,\mu_i)=
\sum^{2N}_{j=1}\int^{\infty}_{0}G(x,\mu_i;x',\mu_j)q(x',\mu_j)\,\mathrm{d}x',
\end{equation}
where the Green's function defined for each $p$ satisfies
\begin{align}
\label{analytical solution3_8}
(\eta+\mu_i)\frac{\partial}{\partial x}G(x,\mu_i;x_0,\mu_{i_0})+
\mu_tG(x,\mu_i;x_0,\mu_{i_0})
&=
\frac{\mu_s}{2}\sum^{2N}_{j=1}w_jG(x,\mu_j;x_0,\mu_{i_0})
\nonumber \\
&+
\delta(x-x_0)\delta_{ii_0},
\end{align}
where $\delta_{ii_0}$ is the Kronecker delta. Here the boundary condition is given by
\begin{equation}
G(0,\mu_i;x_0,\mu_{i_0})=0,\quad\mu_i>\mu_{N_{\eta}}.
\end{equation}

Let us consider the following homogeneous equation to calculate the Green's function with ADO (\citealt{Siewert_1999,Barichello_2000,Barichello_2001}). Note that the equation will be solved for each $p$.
\begin{equation}
\label{analytical solution4_1}
\left((\eta+\mu_i)\frac{\partial}{\partial x}+\mu_t\right)\hat{\psi}(x,\mu_i)=\frac{\mu_s}{2}\sum^{2N}_{j=1}w_j\hat{\psi}(x,\mu_j).
\end{equation}
We note that $\hat{\psi}$ depends on $p$ through $\mu_t$. With separation of variables, we can write $\hat{\psi}$ as
\begin{equation}
\label{analytical solution4_2}
\hat{\psi}(x,\mu_i)=\phi(\nu,\mu_i)\mathrm{e}^{-x/\nu},
\end{equation}
where $\nu$ is the separation constant. The function $\phi(\nu,\mu_i)$ satisfies the normalization condition,
\begin{equation}
\label{analytical solution4_3}
\sum^{2N}_{i=1}w_i\phi(\nu,\mu_i)=
\sum^{N}_{i=1}w_i\left(\phi(\nu,\mu_i)+\phi(\nu,-\mu_i)\right)=1.
\end{equation}
For the integral equation in which the sum on the right-hand side of (\ref{analytical solution4_1}) is replaced by the integral $\frac{\mu_s}{2}\int_{-1}^1\hat{\psi}\,\mathrm{d}\mu$, this $\phi$ is called the singular eigenfunction (\citealt{Case_1960}). However, by the discretization, we obtain
\begin{equation}
\label{analytic_solution_4_5}
\phi(\nu,\mu_i)=\frac{\mu_s\nu}{2}\frac{1}{\mu_t\nu-\mu_i-\eta}.
\end{equation}
assuming $\nu\neq(\mu_i+\eta)/\mu_t$. It is known that $\mu\neq\mu_i/\mu_t$ if $\eta=0$ and $p$ is real (\citealt{Siewert_1999}). We can show the following orthogonality relation.
\begin{equation}
\label{analytic_solution_4_8}
\sum^{2N}_{i=1}w_i(\mu_i+\eta)\phi(\nu,\mu_i)\phi(\nu',\mu_i)=
\mathcal{N}(\nu)\delta_{\nu\nu'},
\end{equation}
where
\begin{equation}
\label{analytic_solution_4_9}
\mathcal{N}(\nu)=
\sum^{2N}_{i=1}w_i(\mu_i+\eta)\phi(\nu,\mu_i)^2.
\end{equation}
We can find $2N$ eigenvalues $\nu=\nu_n$ ($n=1,2,\dots,2N$). Moreover there are $N_{\eta}$ eigenvalues with positive real parts. See Appendix \ref{eigen} for the eigenvalues.

Taking the fact that $G$ vanishes as $x\to\infty$, we can write
\begin{align}
\label{analytic_solution_5_1}
G(x,\mu_i;x_0,\mu_{i_0})
&=
G_0(x,\mu_i;x_0,\mu_{i_0})
\nonumber \\
&+
\sum_{\Re{\nu_n}>0}B(\nu;x_0,\mu_{i_0})\phi(\nu,\mu_i)\mathrm{e}^{-x/\nu},
\end{align}
where coefficients $B(\nu;x_0,\mu_{i_0})$ will be determined later. The free-space Green's function $G_0$ is calculated as
\begin{equation}
\label{analytic_solution_5_6}
G_0(x,\mu_i;x_0,\mu_{i_0})=
\pm\sum_{\pm\Re{\nu_n}>0}\frac{w_{i_0}}{\mathcal{N}(\nu_n)}
\phi(\nu_n,\mu_{i_0})\phi(\nu_n,\mu_i)\mathrm{e}^{-(x-x_0)/\nu_n},
\end{equation}
where upper signs are chosen for $x>x_0$ and lower signs are used for $x<x_0$. Since $G$ vanishes at $x=0$ for $\mu_i>\mu_{N_{\eta}}$, we have from (\ref{analytic_solution_5_1}),
\begin{equation}
\label{analytic_solution_5_7}
\sum_{\Re{\nu_n}>0}B(\nu;x_0,\mu_{i_0})\phi(\nu_n,\mu_i)=
\sum_{\Re{\nu_n}<0}\frac{w_{i_0}}{\mathcal{N}(\nu_n)}
\phi(\nu_n,\mu_{i_0})\phi(\nu_n,\mu_i)\mathrm{e}^{x_0/\nu_n}.
\end{equation}
Let us multiply (\ref{analytic_solution_5_7}) by $\exp(-x_0\mu_t/(\eta+1))$, then integrate both sides with $x_0$ and take the sum with respect to $i$. We obtain
\begin{equation}
\label{analytic_solution_6_1}
\sum_{\Re{\nu_n}>0}E(\nu_n)\phi(\nu_n,\mu_i)=
\sum_{\Re{\nu_n}<0}
\frac{-\nu_n(\eta+1)}{\mathcal{N}(\nu_n)(\eta+1-\nu_n\mu_t)}\phi(\nu_n,\mu_i),
\end{equation}
where
\begin{equation}
\label{analytic_solution_5_9}
E(\nu_n)=\sum_{i_0=1}^{2N}\int_0^{\infty}B(\nu_n;x_0,\mu_{i_0})
\mathrm{e}^{-\mu_tx_0/(\eta+1)}\,\mathrm{d}x_0.
\end{equation}
We can numerically obtain $E(\nu_n)$ from the linear system (\ref{analytic_solution_6_1}). Finally, we obtain
\begin{align}
\label{analytic_solution_6_3}
\hat{\psi_s}(x,\mu_i)
&=
\frac{n_0\mu_s}{2p}\sum_{\Re{\nu_n}>0}\left[
\frac{\nu_n}{\mathcal{N}(\nu_n)}\phi(\nu_n,\mu_i)
\frac{\eta+1}{\eta+1-\nu_n\mu_t}
\left(\mathrm{e}^{-\mu_tx/(\eta+1)}-\mathrm{e}^{-x/\nu_n}\right)
\right]
\nonumber \\
&+
\frac{n_0\mu_s}{2p}\sum_{\Re{\nu_n}<0}\left[
\frac{\nu_n}{\mathcal{N}(\nu_n)}\phi(\nu_n,\mu_i)
\frac{\eta+1}{\eta+1-\nu_n\mu_t}\mathrm{e}^{-\mu_tx/(\eta+1)}
\right]
\nonumber\\
&+
\frac{n_0\mu_s}{2p}\sum_{\Re{\nu_n}>0}E(\nu_n)\phi(\nu_n,\mu_i)
\mathrm{e}^{-x/\nu_n}.
\end{align}
In this way, $\hat{\psi}_s$ is computed using ADO.

The Laplace transform of the particle number density $\hat{n}(x,p)$ is obtained as
\begin{align}
\label{analytic_solution_6_4}
\hat{n}(x,p)&=
\int^1_{-1}\hat{\psi}(x,\mu,p)\,\mathrm{d}\mu
\nonumber\\
&=
\int^1_{-1}\hat{\psi_b}(x,\mu,p)\,\mathrm{d}\mu+
\int^1_{-1}\hat{\psi_s}(x,\mu,p)\,\mathrm{d}\mu
\nonumber\\
&=\frac{n_0}{p}\mathrm{e}^{-x\mu_t/(\eta+1)}
\nonumber\\
&+
\frac{n_0\mu_s}{2p}\sum_{\Re{\nu_n}>0}\left[
\frac{\nu_n}{\mathcal{N}(\nu_n)}\frac{\eta+1}{\eta+1-\nu_n\mu_t}\left(
\mathrm{e}^{-x\mu_t/(\eta+1)}-\mathrm{e}^{-x/\nu_n}\right)\right]
\nonumber\\
&+
\frac{n_0\mu_s}{2p}\sum_{\Re{\nu_n}<0}\left[
\frac{\nu_n}{\mathcal{N}(\nu_n)}\frac{\eta+1}{\eta+1-\nu_n\mu_t}
\mathrm{e}^{-x\mu_t/(\eta+1)}\right]
\nonumber\\
&+
\frac{n_0\mu_s}{2p}\sum_{\Re{\nu_n}>0}E(\nu_n)\mathrm{e}^{-x/\nu_n}.
\end{align}
We note that the particle number density $n(x,t)$ is given by
\begin{equation}
n(x,t)=
\frac{1}{2\pi\mathrm{i}}\int^{\gamma+\mathrm{i}\infty}_{\gamma-\mathrm{i}\infty}\mathrm{e}^{pt}\hat{n}(x,p)\,\mathrm{d}p.
\end{equation}
The Bromwich integral in the inverse Laplace transform is numerically evaluated by the trapezoidal rule. In general, deformation of the contour can be considered for the inversion (\citealt{Weideman_2007}). We found $\gamma=0.04$ is suitable. Hence, the concentration is obtained as
\begin{align}
C(x,t)
&=
\frac{\alpha n_0\mu_s}{4\pi\mathrm{i}}\int^{\gamma+\mathrm{i}\infty}_{\gamma-\mathrm{i}\infty}\frac{\mathrm{e}^{pt}}{p}
\Biggl\{\frac{2}{\mu_s}\mathrm{e}^{-x\mu_t/(\eta+1)}
\nonumber\\
&+
\sum_{\Re{\nu_n}>0}\left[
\frac{\nu_n}{\mathcal{N}(\nu_n)}\frac{\eta+1}{\eta+1-\nu_n\mu_t}\left(
\mathrm{e}^{-x\mu_t/(\eta+1)}-\mathrm{e}^{-x/\nu_n}\right)\right]
\nonumber\\
&+
\sum_{\Re{\nu_n}<0}\left[
\frac{\nu_n}{\mathcal{N}(\nu_n)}\frac{\eta+1}{\eta+1-\nu_n\mu_t}
\mathrm{e}^{-x\mu_t/(\eta+1)}\right]
\nonumber\\
&+
\sum_{\Re{\nu_n}>0}E(\nu_n)\mathrm{e}^{-x/\nu_n}
\Biggr\}\,\mathrm{d}p.
\end{align}

\section{Column experiments}
\label{sec:column}

Column experiments are often used for the study of transport in random media (\citealt{Cotris_2004}). We here present results of a series of tracer breakthrough experiments conducted in a one-dimensional flow field. Fluid with tracer particles is injected by the peristaltic pump (MP-1000, Eyela). The flow rate of injected tracer solution is controlled by the peristaltic pump. Ultra pure water is used. First, water is injected for 24 hours or more until steady flow field is achieved. Next, tracer solution is injected to displace the fresh water. The discharged solution is collected by the fraction collector (CHF161RA, Advantec) at regular intervals. Column experiments are performed at room temperature of $25$--$26\,{}^{\circ}{\rm C}$.

Tracer experiments are conducted until the discharge concentration becomes equal to the influent concentration. For the sake of the accuracy of measurements, we correct the measurement time error from the tube volume and inlet volume of the column by subtracting the time lag due to the switchover from the breakthrough time.

\subsection{Column experiments with non-adsorbed solute}

First we use non-adsorbed solute. We use two different sizes of glass beads made with Soda glass; (Run 1) smaller beads with diameter ranging from $0.177\,{\rm mm}$ to $0.250\,{\rm mm}$ and (Run 2) larger beads with diameter from $0.500\,{\rm mm}$ to $0.710\,{\rm mm}$. The column length is $20.0\,{\rm cm}$ with the internal diameter of the column $2.9\,{\rm cm}$. In order to prevent entrainment of air, we carry out beads packing under saturated conditions, in which the water surface is kept at a constant height from the top of the beads layer. We stirrer beads while filling in order to remove small bubbles attached to the beads.

The inlet and the outlet ends of the column are separated from the porous medium of a glass filter with the pore size $10\,{\rm \mu m}$, originally for the purpose of the liquid chromatography. The glass filter serves to enhance tracer mixing and to dampen small flow pulses on injection. The tracer concentration is measured by the EC meter (LAQUAtwin-EC-33B, Horiba). The precision of the meter is $\pm 2\%$. The temperature compensation circuit is installed.

We prepare the NaCl solution of $0.1\%$ with the salt of purity $99.5\%$ (Wako). Injection of the solution is done by the peristaltic pump. The measured EC values are given by the unit of Sv/cm, we made the calibration curve and obtain the concentration.

\subsubsection{Run 1: small beads}
\label{exp_run1}

The column is filled with small beads. The measured porosity is $\eta=0.386$. The bed height is $18.0\,{\rm cm}$ with section area $6.61\,{\rm cm}^2$. The flow rate is $0.958\,{\rm cm}^3/{\rm min}$.

Figure \ref{fig:run1} shows the chloride breakthrough curve from Run 1. The vertical axis shows the relative concentration $C/C_0$ and the horizontal axis is time $t$. The measured values from the column experiment (black dots) are compared to $C/C_0$ calculated by the linear Boltzmann equation (LBE) (red solid line) and by the ADE model (\citealt{Ogata_1961}) (blue dashed line). The parameters in both equations are determined by nonlinear least-squares fitting (i.e., the minimization of the sum of squares of the residuals between the experimental and theoretical curves); they are obtained as $\sigma_a=1.0\times 10^{-8}\,{\rm min}^{-1}$, $\sigma_s=5.1645\,{\rm min}^{-1}$, $v_0=5.3073\,{\rm cm}/{\rm min}$, $u=1.6445\,{\rm cm}/{\rm min}$, $\beta=0.09130$ for LBE, and $u=1.2886\,{\rm cm}/{\rm min}$, $D=1.8379\,{\rm cm}^2/{\rm min}$ for ADE. In Fig.~\ref{fig:run1}, experimental results agree well with both LBE and ADE. Figure \ref{fig:run1_log} shows semi-log plots of $1-C/C_0$ as a function of $t$.

\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]{run1_nacl.eps}
\caption{
(Color online) The relative concentration $C/C_0$ for Run 1 (black dots) is plotted with calculated values from LBE (red solid line) and ADE (blue dashed line).
}
\label{fig:run1}
\end{figure}

\begin{figure}
\centering
\includegraphics[clip,width=0.8\textwidth]{run1log_nacl.eps}

\caption{
(Color online) Same as Fig.~\ref{fig:run1} but $1-C/C_0$ is plotted in the semi-log graph.
}
\label{fig:run1_log}
\end{figure}

\subsubsection{Run 2: large beads}
\label{exp_run2}

The column is filled with large beads. The measured porosity is $\eta=0.367$. (This value of the porosity is smaller than that of Run 1. This may be due to small bubbles attached on the beads surface during Run 1.) The bed height is $18.5\,{\rm cm}$ with section area $6.61\,{\rm cm}^2$. The flow rate is $0.875\,{\rm cm}^3/{\rm min}$.

Figure \ref{fig:run2} for Run 2 is the same as Fig.~\ref{fig:run1} but large glass beads are used. By nonlinear least-squares fitting, the parameters are obtained as $\sigma_a=1.0\times 10^{-8}\,{\rm min}^{-1}$, $\sigma_s=4.6778\,{\rm min}^{-1}$, $v_0=8.1753\,{\rm cm}/{\rm min}$, $u=2.4999\,{\rm cm}/{\rm min}$, $\beta=0.08579$ for LBE, and $u=1.6189\,{\rm cm}/{\rm min}$, $D=3.8519\,{\rm cm}^2/{\rm min}$ for ADE. In Fig.~\ref{fig:run2}, experimental results are described well by LBE and ADE. The long-time behavior is shown in Fig.~\ref{fig:run2_log}, in which $1-C/C_0$ is plotted.

\begin{figure}
\centering
\includegraphics[clip,width=0.8\textwidth]{run2_nacl.eps}
\caption{
(Color online) The relative concentration $C/C_0$ for Run 2 (black dots) is plotted with calculated values from LBE (red solid line), and ADE (blue dashed line).
}
\label{fig:run2}
\end{figure}

\clearpage
\begin{figure}
\centering
\includegraphics[clip,width=0.8\textwidth]{run2log_nacl.eps}
\caption{
(Color online) Same as Fig.~\ref{fig:run2} but $1-C/C_0$ is plotted in the semi-log graph.
}
\label{fig:run2_log}
\end{figure}

\subsection{Column experiments with adsorbed solute}

Here we use adsorptive solute (zinc solution). The filling material is the standard sand (Tohoku silica sand No. 4, Kitanihon Sangyo). The median diameter of the sand is 750$\mu$m. As a preparation, we eliminate the organic matter that may have been contained in the sand by soaking it in ${\rm HNO}_3$ solution. The zinc solution is 2 ppm. We set a filter on the top of the sand bed made with glass wool. We also put the same filter at the bottom of the column.

A short column of length $12.0\,{\rm cm}$ with internal diameter $3.1\,{\rm cm}$ was used.
We perform a blank test beforehand to make sure that the zinc is not absorbed on the surface of the column wall. The concentration is measured with the atomic absorption photometer (Z-2300, Hitachi High-Technologies), the compressor (SC820, Koki Holdings), and the neo cool circulator (CF700, Yamato).

\subsubsection{Run 3}
\label{exp_run3}

The measured porosity is $0.295$ and the flow rate is $9.34\,{\rm cm}^3/{\rm h}$. The bed height is $10.7\,{\rm cm}$ and the section area is $7.54\,{\rm cm}^2$.

Figure \ref{fig:run3} shows the relative concentration $C/C_0$ of the zinc breakthrough curve for Run 3 (black dots) with calculated values from the linear Boltzmann equation (LBE) (red solid line) and the Ogata-Banks ADE model (blue dashed line). By nonlinear least-squares fitting, we obtain $\sigma_a=1.0\times 10^{-7}\,{\rm h}^{-1}$, $\sigma_s=2.8134\,{\rm h}^{-1}$, $v_0=5.0663\,{\rm cm}/{\rm h}$, $u=1.9876\,{\rm cm}/{\rm h}$, $\beta=0.1739$ for LBE, and $u=1.1406\,{\rm cm}/{\rm h}$, $D=4.3864\,{\rm cm}^2/{\rm h}$. The dashed curve for ADE partially deviates from the measured values. The semi-log plot in Fig.~\ref{fig:run3_log} shows a clear discrepancy between the experimental result and ADE.

\begin{figure}
\centering
\includegraphics[clip,width=0.8\textwidth]{run3_zn.eps}
\caption{
The relative concentration $C/C_0$ for Run 3 (black dots) is plotted with calculated values by LBE (red solid line) and by ADE (blue dashed line).
}
\label{fig:run3}
\end{figure}

\begin{figure}
\centering
\includegraphics[clip,width=0.8\textwidth]{run3log_zn.eps}
\caption{
Same as Fig.~\ref{fig:run3} but $1-C/C_0$ is plotted in the semi-log graph.
}
\label{fig:run3_log}
\end{figure}

\subsubsection{Run 4}
\label{exp_run4}

Figure \ref{fig:run4} for Run 4 is the same as Fig.~\ref{fig:run3} but we repeat another run for the standard sand and adsorbed solute. The measured porosity is $0.289$ and the flow rate is $11.25\,{\rm cm}^3/{\rm h}$. The bed height is $10.7\,{\rm cm}$ and the section area is $7.54\,{\rm cm}^2$. By nonlinear least-squares fitting, we obtain $\sigma_a=1.0\times 10^{-7}\,{\rm h}^{-1}$, $\sigma_s=2.6773\,{\rm h}^{-1}$, $v_0=4.1166\,{\rm cm}/{\rm h}$, $u=1.7999\,{\rm cm}/{\rm h}$, $\beta=0.2091$ for LBE, and $u=1.1384\,{\rm cm}/{\rm h}$, $D=3.9814\,{\rm cm}^2/{\rm h}$ for ADE. The discrepancy of ADE is even more apparent. The semi-log plot in Fig.~\ref{fig:run4_log} shows a discrepancy between the experimental result and ADE.

\begin{figure}
\centering
\includegraphics[clip,width=0.8\textwidth]{run4_zn.eps}
\caption{
The relative concentration $C/C_0$ for Run 4 (black dots) is plotted with calculated values by LBE (red solid line) and by ADE (blue dashed line).
}
\label{fig:run4}
\end{figure}

\begin{figure}
\centering
\includegraphics[clip,width=0.8\textwidth]{run4log_zn.eps}
\caption{
Same as Fig.~\ref{fig:run4} but $1-C/C_0$ is plotted in the semi-log graph.
}
\label{fig:run4_log}
\end{figure}

\clearpage
\section{Conclusions}
\label{sec:concl}

Not only for field observation reported in studies such as Adams and Gelhar (\citealt{Adams_1992}) and Benson and his collaborators (\citealt{Benson_2001,Benson_2000}), anomalous diffusion also appears in lab-scale experiments (\citealt{Cotris_2004}). In Sec.~\ref{sec:column}, ADE does not precisely reproduce the breakthrough curves even for beads, which form a medium material with non-adsorbed solution. Williams proposed to use the linear Boltzmann equation for transport in porous media (\citealt{Williams_1992a,Williams_1992b,Williams_1993a,Williams_1993b}). It was shown that the proposed model of the linear Boltzmann equation reproduces the whole breakthrough curves of experimental results including the non-monotonic decay of $1-C/C_0$ in the tail.

The diffusion approximation holds when the propagation distance of tracer particles is sufficiently larger than the transport mean free path $\ell^*$, which is given in (\ref{LBE_ADE_9}). In the case of light transport, the photon diffusion equation starts to work when the propagation distance becomes $10$ times larger than the transport mean free path (\citealt{Yoo_1990}). Using the relation (\ref{LBE_ADE_8}), we can calculate the diffusion coefficient, which is denoted by $D'$, from parameters in the linear Boltzmann equation. We let $D$ denote the fitted diffusion coefficient in the advection-diffusion equation. Let us compare Run 1 and Run 2. For Run 1 in Sec.~\ref{exp_run1}, $\ell^*=1.03$ and $D'=1.818$, whereas $D=1.8379$. The relative difference is $|D-D'|/D=0.0108$. For Run 2 in Sec.~\ref{exp_run2}, $\ell^*=1.75$ and $D'=4.762$, whereas $D=3.8519$. The relative difference is $|D-D'|/D=0.236$. In either case, the propagation distance of $18\,{\rm cm}$ or $18.5\,{\rm cm}$ is much larger than $\ell^*$ and the transport is expected to be in the diffusion regime. Indeed, Fig.~\ref{fig:run1} and Fig.~\ref{fig:run2} show that both LBE and ADE reproduce the experimentally obtained breakthrough curves. Next, let us also look at Run 3 and Run 4. For Run 3 in Sec.~\ref{exp_run3}, $\ell^*=1.80$ and $D'=3.041$, whereas $D=4.3864$. The relative difference is $|D-D'|/D=0.307$. For Run 4 in Sec.~\ref{exp_run4}, $\ell^*=1.54$ and $D'=2.110$, whereas $D=3.9814$. The relative difference is $|D-D'|/D=0.470$. The propagation distance of $10.7\,{\rm cm}$ is not as long as the distance in Run 1 and Run 2 in terms of the transport mean free path. The discrepancy of the curve for ADE in Fig.~\ref{fig:run3} and Fig.~\ref{fig:run4} might also attribute to anisotropic scattering for the sand. The transport mean free path becomes larger when the scattering in LBE becomes anisotropic. Moreover the discrepancy may come from the fact that the surface of the sand is rough and the particle speed cannot be modeled by a constant $v_0$.

Although we assumed isotropic scattering in (\ref{proposed_model}), it is possible to take anisotropic scattering into account and write the integral term in a more general form as $\sigma_s\int_{-1}^1p(\mu,\mu')\psi(x,\mu',t)\,\mathrm{d}\mu'$, where $p(\mu,\mu')$ is the probability that tracer particles moving in the direction $\mu'$ change their direction to the direction $\mu$ by scattering. We can similarly develop the analytical discrete ordinates method for anisotropic scattering by writing $p(\mu,\mu')=\frac{1}{2}\sum_l\beta_lP_l(\mu)P_l(\mu')$ with some constant coefficients $\beta_l$ and Legendre polynomials $P_l$ ($l=0,1,\dots$) (\citealt{Barichello_2011,Siewert_2000}).

In this paper we have shown that the mass transport in column experiments, which show anomalous diffusion, obeys the linear Boltzmann equation. Furthermore the linear Boltzmann equation in the transport regime appears at the mesoscopic scale and the advection-diffusion equation in the diffusion regime is derived from the linear Boltzmann equation in the asymptotic limit when the propagation distance is sufficiently larger than $\ell^*$ at the macroscopic scale. The physical origin of anomalous diffusion in a heterogeneous random medium with more complex structure of fissures is still an open problem.

\begin{acknowledgements}
The seed of this research was planted on the occasion of the Study Group Workshop (Department of Mathematical Sciences, The University of Tokyo, December 2010), which is greatly appreciated. The research was restarted with the support of the Focusing Collaborative Research by the Interdisciplinary Project on Environmental Transfer of Radionuclides (University of Tsukuba and Hirosaki University). The research was partially supported by the JSPS A3 foresight program: Modeling and Computation of Applied Inverse Problems. MM acknowledges support from Grant-in-Aid for Scientific Research (17H02081,17K05572) of the Japan Society for the Promotion of Science. YH acknowledges support from Grant-in-Aid for Scientific Research (15H05740,17H01478) of the Japan Society for the Promotion of Science.
\end{acknowledgements}