\section{Introduction.}

The spin-down of an isolated neutron star (NS) on secular timescales is mainly caused by rotational energy losses due to electromagnetic radiation, pulsar winds, or
gravitational radiation.
A measurable quantity closely related to the rotational evolution of pulsars is the braking index $n$,
defined by assuming that the star spins down according to a power law $\dot{\Omega}= -K \Omega^n$,
where $\Omega$ is the  NS angular velocity.
For each of these dominant rotational energy loss mechanisms, the proportionality
constant $K$ hides different dependences on the star radius, moment of inertia, magnetic field strength,
and angle between rotation and magnetic axis.
If all these quantities are constant in time, the magneto-dipole spindown mechanism
predicts a braking index $n=3$, but variations in time of any of these quantities may cause
departures from this canonical value.
Unfortunately, the accurate determination of the second derivative of the frequency, needed to estimate the  braking
index, is not always possible because it is affected by glitches and other short-term timing irregularities.

At present, eight pulsars have sufficiently steady rotations that stable values of their braking index are generally accepted \citep{Lyne93,Lyne96,Mid2006,Liv2007,Liv2011,Welt2011,Espinoza11}, among which the most recent case is
PSR J1734-3333, which has $n=0.9\pm0.2$ significantly below 3 \citep{Espinoza11}.
All these cases show the same trend: they are all young pulsars (Vela is the oldest among these pulsars
with $t\approx 12$ kyr) and they all have $n<3$.
We note that gravitational wave emission predicts $n=5$, but
it is only efficient during the first minutes or hours of a NS life, when rotation
is sufficiently fast and the mass quadrupole moment large enough (see e.g., \cite{HJA06} and references therein).
It certainly does not contribute significantly to the timing noise of pulsars older than $100$ yr.

Selecting a sample of 127 pulsars from the ATNF Pulsar Catalogue \citep{ATNF},
for which the quoted errors in  the second derivative of the spin frequency ($\ddot{\nu}$) are less than ten per cent,
\cite{Urama2006}
found a strong correlation of $\ddot{\nu}$ with $\dot{\nu}$, independent of the sign of $\ddot{\nu}$. They suggested
that this trend can be accounted for by small stochastic deviations in the spin-down torque that are directly
proportional (in magnitude) to the spin-down torque itself.
Another point discussed in the literature is that some of the old pulsars ($>10^6$ yr) have braking indices with absolute values exceeding $|n|=10^4$.
The occurrence of very high braking indices of both signs has been considered in the context of internal frictional instabilities occurring between the crust and the superfluid, almost independently of the evolution of the neutron star magnetic field \citep{SM95}. However, this applies only for old neutron stars ($\tau \gtrsim 2 \times 10^7$ yr) and appears as extremely short term events oscillating about the canonical value $n=3$.
\cite{BT10} proposed another explanation for the observed distribution with very high positive and negative braking
indices by studying the effect of nondipolar magnetic field components and neutron star precession
on magnetospheric electric current losses.
These large $n$ should be observable over relatively long periods of $10^3$-$10^4$ yr.
Another possibility that can explain the observed variability of braking indices is
the time-evolution of conductivity in the magnetosphere \citep{Spit2012}, which also has
implications for the spin-down of intermittent pulsars and subpulse drift phenomena \citep{Lyne2009}.

In this paper we focus on the imprint that the time-evolution of the internal magnetic field has
on the timing noise. Our goal is to estimate the contribution of this variability (expected in realistic evolutionary models)
to the braking index at different epochs. In particular, we discuss the possible reemergence of a magnetic field
initially submerged by hypercritical accretion during the supernova explosion, and the evolution of the crustal
magnetic field under the combined action of Ohmic dissipation and Hall drift during the first few million years of a NS life.

\section{Basic equations.}

The spin--down behavior of a rotating neutron star is governed by the energy balance equation
relating the loss of rotational energy,  ${E}_{\rm{rot}} = I\Omega^2/2$,
where $I$ is the moment of inertia of the neutron star and $\Omega$ the angular velocity, to the
energy loss rate $ \dot{E}$, given by magneto-dipole radiation, wind, gravitational radiation, or others:
\begin{equation}
\dot{E}_{\rm{rot}} = I \Omega \dot{\Omega} \approx \dot{E}.
\hspace*{-5mm}\label{equ:spin-down}
\end{equation}

The standard way to define the braking index is
\begin{equation}
 n = \frac{\ddot{\Omega}\Omega}{\dot{\Omega}^2} = \frac{\ddot{\nu} \nu}{\dot{\nu}^2} =
 2 - \frac{\ddot{P}P}{\dot{P}^2} \;,
\hspace*{-5mm}\label{equ:BIstandard}
\end{equation}
where  $\nu$ is the spin frequency, and we denote by $P=1/\nu$ the period.
Under the usual assumption of a power-law rotational evolution
\begin{equation}
\dot{\Omega} = - K \Omega^{n}\,,
\label{equ:omegadot-hat}
\end{equation}
the power-law index coincides with the braking index.

For the particular case of a rotating dipole in vacuum, the well--known Larmor formula returns
\begin{equation}
\dot{E}= \frac{B_0^2R^6\Omega^4}{3c^3} \sin{\alpha}^2\;,
\label{equ:mdr}
\end{equation}
where $B_0$ is the dipolar component of the surface magnetic field at the pole,  $R$ denotes the neutron star radius,
$\alpha$ is the angle between the rotational and the magnetic axis, and $c$ is the speed of light.
We neglect the contribution to the spindown of higher order multipoles,
since they have significantly shorter ``lever arms" than the dipolar mode.

Combining Eqs. (\ref{equ:spin-down}) and (\ref{equ:mdr}), one obtains
\begin{equation}
\dot{\Omega} = - K \Omega^3
\label{equ:omegadot}\;,
\end{equation}
where $K=f_\alpha B_0^2R^6/3Ic^3$.  For simplicity, but without loss of generality, we have omitted the
geometric factor that carries the dependence on the inclination angle, and we
assume hereafter that it is a correction $f_\alpha$ factorized in $K$, which carries information about the
particular physical process that governs energy losses.
The functional dependence of $K$ for all other magnetic processes is the same $\propto B_0^2\Omega^4R^6$.
Differences in the radiation mechanism are included in the factor $f_\alpha$. While
magnetospheric current losses scale as $\cos^2\alpha$ \citep{BN07},
magneto-dipole radiation losses scale as $\sin^2\alpha$ for vacuum or
$\frac{3}{2}(1+\sin^2\alpha)$ for force-free magnetospheres \citep{Spitko2006}.
The most recent resistive solutions for pulsar magnetospheres \citep{Spit2012} fit the spin-down luminosity with a prefactor
of the order of unity that also depends on the maximum potential drop along field lines in the corotating frame.

Eq.~(\ref{equ:omegadot}) can also be cast in the usual form
\begin{equation}
P \dot{P} = {\cal{K}} B_0^2(t)
\label{equ:spin-down-B(t)}
\end{equation}
where for a standard neutron star ($R=10^6$cm, $I=10^{45}$gcm$^2$) the constant
${\cal{K}} \approx 10^{-39}$ cm s$^3$g$^{-1}$. Here we have explicitly written $B_0(t)$ to show that the
magnetic field is the only quantity that we allow to vary with time.

Now we consider what happens when a secular variation of the magnetic field is allowed. Deriving $\ddot{P}$ from
Eq.~(\ref{equ:spin-down-B(t)}), we have
\begin{equation}
\ddot{P}={\cal{K}}BP^{-2}\left(2\dot{B_0}P\,\,-\,\,B_0\dot{P}\right)~,
\end{equation}
and the braking index in a general case of a time-dependent magnetic field can be simply expressed as follows
\begin{equation}
n=3-4 \frac{\dot{B_0}}{B_0} \tau_c \equiv 3-4 \frac{\tau_c}{\tau_B}~,
\hspace*{-5mm}\label{equ:hatn}
\end{equation}
where $\tau_c = P/2 \dot P$ is the characteristic age, and we defined
the magnetic field evolution timescale
$\tau_B \equiv B_0/\dot{B_0}$.
This last equation shows that any
variation of the dipolar surface magnetic field strength results in a deviation from the $n=3$
standard value,  which is obviously recovered for a constant magnetic field ($\dot{B_0} = 0$).
For an increasing $B_0$ we will always obtain ${n} < 3$, while ${n} > 3$ is the signature
of a decreasing $B_0$.

\section{Magnetic field evolution scenarios.}

The precise calculation of the time variation of $B_0$ requires solving the problem of the coupled magneto-thermal
evolution of a neutron star \citep{PMG2009,VPM2012}.
However, before showing results from simulations, one can qualitatively analyze
of the most important physical process at different ages.
We now discuss three possible scenarios.

\subsection{Amplification or reemergence of the dipolar surface magnetic field.}

The generic $n< 3$ observed without exception for young pulsars can be
caused by the rediffusion of the magnetic field submerged in the crust during the supernova fallback episode
\citep{young95,muslimov95,GPZ99}. The submergence of the magnetic field has also been studied
in recent MHD simulations, with interesting implications for gravitational wave emission \citep{Vigelius2009,Vigelius2010}.
The rediffusion timescale depends essentially on the total amount of accreted matter.
For total accreted masses between $10^{-4}$-$10^{-3} M_\odot$, screening currents are dissipated on
$\tau_{\rm rediff} \sim 10^3$-$10^4$ yr, as recent calculations have confirmed \citep{Ho2011,Vigano12}.
The same idea has been applied to the reemergence of magnetic fields in accreting binary
systems with similar results \citep{CZB2001}, and evidence for
fast magnetic field evolution in an accreting millisecond pulsar has been reported \citep{Patruno2012}.

Alternative mechanisms are time variations in the angle between the rotational and magnetic axes
\citep{LE1997,RZC1998},
or the thermoelectric field generation that may proceed in the crust and envelope
of young pulsars if a sufficiently strong temperature gradient is present  \citep{ULY86,WG96}.
This process  is limited by the condition that the surface temperature of the neutron star should not be lower than
$3\times 10^6$ K \citep{WG96}, which corresponds to about $1000$ yr in the standard cooling scenario.
In our simulations we included the reemergence of a screened magnetic field, but our present version of the
code does not include magnetic field generation by thermoelectric effect, or other possible mechanisms
such as magnetic flux expulsion from the superconducting core.

\subsection{Ohmic decay.}

Ohmic dissipation of the magnetic field has been thoroughly investigated in the past
by many groups \citep{HUY1990,PGZ00,Tau2001}.
The usual definition of the Ohmic decay time is
\begin{equation}
\tau_{\rm{Ohm}}= \frac{4\pi\sigma L^2}{c^2}\;,
\hspace*{-5mm}\label{equ:tauohm}
\end{equation}
where $L$ denotes the scale-length of magnetic field variations and $\sigma$ the electric conductivity.
In a neutron star crust, $\sigma$ is dominated by electronic transport and depends on the
electron density, the crustal temperature, and the impurity concentration within the crust.
Since the electron density varies over about four orders of magnitude in the crust and the temperature decreases by about two to three orders of magnitude in a pulsar's lifetime, the electric conductivity may vary both in space and time by many orders of magnitude (see Fig. 1 in \cite{PG07}). Therefore, to assume a uniform Ohmic decay time $\tau_{\rm{Ohm}}$, independent of the location of currents and the pulsar age is, generally speaking, quite misleading:
Ohmic decay in a neutron star crust cannot be described by a single exponential law
(see e.g. \cite{PGZ00} for a qualitative discussion).
Since the conductivity increases with decreasing temperature,  Ohmic diffusion becomes
an increasingly slow process in older pulsars, beyond the transition from neutrino dominated
to photon cooling at $\gtrsim 10^5$ yr. The only general fact that we can expect is that,
for low-field middle-aged neutron stars ($10^4$-$10^5$ yr),
this is expected to be the dominant process, thus resulting in $n>3$.

\subsection{Hall drift oscillatory modes.}
The Hall drift, through its nonlinear dependence on the magnetic field,
has an influence on the magnetic field evolution either for magnetar conditions ($B>10^{14}$ G), or alternatively when the
electrical resistivity becomes very low.
Under certain circumstances the Hall drift can drain magnetic energy out of the dipolar mode and redistribute it into smaller scale ones. This sometimes causes an apparently rapid decrease of $B_0$. But the
opposite effect may also occur, i.e., pumping energy from an internal strong toroidal field to the dipolar poloidal component.
Both scenarios would happen on the Hall timescale
\begin{equation}
\tau_{\rm{Hall}}= \frac{4\pi n_e e L^2}{cB}\;,
\hspace*{-5mm}\label{equ:tauHall}
\end{equation}
where $n_e$ is the electron density and $e$ the elementary charge. For typical values in NSs we have
$\tau_{\rm Hall} \sim (10^4$-$10^6) \frac{10^{14}{\rm G}}{B_0}$ yr.

In situations of quasi-equilibrium, oscillatory modes with magnetic energy exchange between the
different crustal field modes are expected. The temporal evolution of the polar surface magnetic field can then be
approximated by
\begin{equation}
B_0(t)= B_0 + \delta B\sin \left( \frac{2\pi t}{\tau_{\rm Hall}}\right)\;,
\hspace*{-5mm}\label{equ:tauHall2}
\end{equation}
and from Eq.~(\ref{equ:hatn}) we have
\begin{equation}
n=3 -4\frac{\dot{B_0}}{B_0} \tau_c \approx 3-8\pi \frac{\tau_c}{\tau_{\rm{Hall}}} \frac{\delta B}{B_0} \cos\left( \frac{2\pi t}{\tau_{\rm Hall}}\right)\;.
\hspace*{-5mm}\label{equ:tauHall3}
\end{equation}

These oscillations are expected in young magnetars with $B_0>10^{14}$ G. It causes
corrections to the canonical braking index $n=3$ of either positive or negative sign, according to Eq.~(\ref{equ:tauHall3}).
These corrections are
expected to be small because of their short $\tau_c$, and to be increasingly
important for objects with smaller $\dot{P}$ (old characteristic ages).
In addition, as the star cools down, the drop in electrical resistivity may activate the Hall term
even for normal pulsars $B_0 \gtrsim 10^{12}$ and the occurrence of oscillatory modes
during a {\it second Hall stage} at late times is a natural outcome.
First estimates suggest that $\delta B/B_0$ could be as large as $\sim 10^{-1}$.
When $\tau_c \gg \tau_{\rm Hall}$,
the second term in Eq.~(\ref{equ:tauHall3}) dominates and the magnetic field oscillatory modes should
result in equally probable positive and negative braking indices with high absolute values. However,
this needs to be confirmed by long term realistic simulations.

\begin{figure}[th]
\centering
\includegraphics[width=0.4\textwidth]{fig1.ps}
\caption{
Top: braking indices (absolute value) as a function of $\tau_c$ for our sample of 118 pulsars.
We represent pulsars with $n<3$ with crosses and pulsars with $n>3$ with diamonds. Note that
all objects marked with crosses and $\tau_c>10^5$ yr have actually $n<0$.
The 8 youngest objects discussed in \cite{Espinoza11} are marked with squares.
Middle:  magnetic field evolution timescale, $\tau_B=4 \tau_c/(3-n)$ for pulsars with $n<3$.
Bottom: Absolute value of $\tau_B$ (they are all negative) for pulsars with $n>3$.  }
\label{fig1}
\end{figure}

\section{Discussion.}

\subsection{Magnetic field evolution in the pulsar population?}

From the whole population in the ATNF Pulsar Catalogue \citep{ATNF}, we extracted a sample of
pulsars for which $\ddot{\nu}$ has a quoted error smaller than 10\%. From this preselected sample,
we excluded all pulsars in binaries and those with very short periods ($P < 15$ ms), likely to be recycled,
and we doubled checked our list with the more recent review by \cite{Hobbs10}, the first large-scale analysis
of pulsar timing noise over timescales $>10$ yr, which led us to remove several more pulsars whose revised values were inconsistent
with the ATNF data, or with larger errors. After this selection, our sample contains 118 radio-pulsars,
about half of which have negative braking indices.

We show in the top panel of Fig.~\ref{fig1} the observed distribution of braking indices as a function of $\tau_c$. We  considered the characteristic magnetic field evolution timescale, $\tau_B= 4 \tau_c/(3-n)$ (see Eq.~(\ref{equ:hatn})) and
separated the sample into two groups, those with positive (middle panel) and negative (bottom panel) $\tau_B$. The eight pulsars discussed in \cite{Espinoza11} are marked with squares.
With all due cautions regarding the uncertainties associated to these measurements, it is worth mentioning some interesting trends visible in this plot:

\begin{itemize}
\item
All young objects seem to have always $n<3$, which can be a hint of an increasing dipolar magnetic field.
\footnote{We use the values collected in \cite{Espinoza11} for this plot.}
\item
All middle age objects ($10^4$-$10^5$ yr), except Vela, that survived to our conservative selection
criteria have negative time derivatives of $B_0$.\footnote{We did not include seven more objects,
that fall in this region, with a quoted error of $\ddot{\nu}$
smaller than 10\% in the ATNF pulsar database. They are not considered in \cite{Hobbs10}, probably because they did not have 10 yr of accumulated data, although their values have not been reported to change.}
\item
For old objects, there is no correlation at all, and there are similar numbers of objects with positive and negative
derivatives of the field. The typical evolution timescales are in the range $10^3$-$10^5$ yr.
\end{itemize}

The strong correlation between $|n|$ and $\tau_c$ seen in the old objects of the top panel simply reflectis the definition of $n$, Eq.~(\ref{equ:hatn}), with the additional piece of information that $\tau_B$ does not seem to be
correlated with $\tau_c$ (see middle and bottom panels).

Note that our criteria automatically select the objects with a dominant contribution of $\ddot{\nu}$ in the timing phase residuals (i.e., the cleanest cubic lines in Fig.~3 of \cite{Hobbs10}, e.g. B0114+58). This introduces a bias toward objects with high braking indices, and against pulsars with residuals dominated by higher order terms (e.g. B0136+57) or quasi-periodic terms (e.g. B1642-03 or B1826-17) in the time-dependent phase \citep{Lyne2010}. In the latter cases, the estimates of $\ddot{\nu}$ are subject to larger uncertainties, and $n$ strongly depends on the analyzed time interval \citep{Hobbs10}. However, we have checked that the general trends do not change if the sample is enlarged by including pulsars with quoted errors in $\ddot{\nu}$ up to 80\%. This simply increases the statistics ($>300$ sources) and includes some objects with lower value of $|n-3|$ (shorter $|\tau_B|$).

\subsection{Expected evolution from theoretical models.}

We performed a series of numerical simulations with the magneto-thermal evolution code presented in \cite{Aguilera2008b,PMG2009}, and \cite{VPM2012}, which includes all relevant microphysical processes and follows the complex feedback between the physical mechanisms presented in a simplified way in the previous section. In addition to Ohmic dissipation and Hall drift effects, which are consistently included in the simulations, we assumed that all models suffer an episode of hypercritical accretion during their formation.

We plot in Fig.~\ref{fig2} the evolutionary tracks of some theoretical models in the $P$-$\dot{P}$ diagram, compared to the pulsars of our sample.
We show results for three initial magnetic field configurations. Two of them are purely dipolar, although during the evolution other multipoles and toroidal field are naturally created. They differ by the magnetic field strengths:
$B_0=10^{13}$ G (model A, dashed lines) and $10^{14}$ G (model B, dotted lines). In each case (models A and B)
we compare results with three values of the total accreted mass.
The third initial configuration (model C, solid line) has an initial dipolar field of $B_0=10^{13}$ G, as model A,
but with an additional octupolar component of strength $3.5\times10^{14}$ G (at the pole).
This model serves to make explicit the effect of the complex initial geometries.
Assuming vacuum external boundary conditions, the
toroidal field has to vanish at the star surface, but the dipolar component
is still coupled to higher order multipoles, which can also interchange energy with the dipolar mode.

\begin{figure}[th]
\centering
\includegraphics[width=0.45\textwidth]{fig2.eps}
\caption{Evolutionary tracks in the $P$-$\dot{P}$ diagram during the first 3 Myr for: model A (dashed) and model B (dotted), both for different values of
accreted mass $M_a=[0.1,5,10] \times 10^{-4} M_\odot$; model C (solid) without accretion.}
\label{fig2}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth]{fig3.eps}
\caption{Braking index ($n$) as a function of age for model A (dashed line), model B with $M_a = 10^{-3} M_\odot$ (thin
solid line), and model C (thick solid line).}
\label{fig3}
\end{figure}

For the eight youngest pulsars, we indicate their predicted movement for the next 2.5 kyr with arrows, assuming the present value of $n$, from \cite{Espinoza11}, remains constant. The direction of the tangent vector to a given track is related to the braking index. For reference, we indicate in the legend on the upper right corner that $n=3$ and $n=1$ imply constant inferred $B_0$, and constant characteristic age, respectively.

The initial period $P_{in}$, which was assumed to be 0.01 s in all cases, only affects the early stage, while $\dot{P} \lesssim P_{in}/t$ (with $t$ being the real age), and $P\simeq P_{in}$. Models with deep submergence of the field
but with different initial periods also have vertical trajectories shifted to the left/right depending on $P_{in}$, and
quickly cross the range of $\dot{P}$ where the bulk of pulsars lie.
At late times, tracks coming from the same model but with different $P_{in}$ are indistinguishable,
typically converging after $t \sim 30 (\frac{P_{in}}{0.01 s}\frac{10^{13}G}{B_{in}})^2$ yr, with $B_{in}$ being
the initial magnetic field.
Therefore, $P_{in}$ has an appreciable long-term effect only in the deep submergence case ($M_a \gtrsim 10^{-3} M_\odot$), for which $B_{in}$ is strongly reduced: in these models, the first few $10^4$ yr are spent in the vertical trajectories,
with $n\ll 0$. Note that in this reemergence phase there is no correlation between the real and characteristic ages.

When reemergence of the field has almost been completed, the trajectories reach the high-$\dot{P}$ region (i.e., largest $B_0$) and progressively bend. The more extreme braking indices of PSR J1734-3333 ($n=0.9\pm0.2$, \cite{Espinoza11}) and PSR J0537-6910 ($n\sim -1.5$, \cite{Mid2006}) would be consistent with the last stage of the reemergence after a deep submergence into the inner crust.
On the other hand, in the shallow submergence models ($M_a \sim 10^{-5}$-$10^{-4} M_\odot$),
after accretion stage $B_0\simeq B_{in}$. These tracks initially run almost along the iso-magnetic lines: pulsars with $n$ slightly less than 3 are compatible with this scenario.

Independently of the early reemergence phase (if any), tracks with the same $B_{in}$ converge at middle-age, and
have slopes corresponding to $n>3$,  characteristic of the slow Ohmic dissipation. During this phase, there is a correlation between real and characteristic ages, with typically $\tau_c$ being a factor of few longer than the real age. Some tracks show visible oscillations produced by the Hall activity when the star is cold enough ($t\gtrsim 10^5$ yr); in particular, model C (solid line) clearly shows that any complex initial geometry may have a distinct
signature on the timing properties of pulsars.

We emphasize that we did not attempt to fit individual objects: our purpose with this sample of models is simply to show that, with reasonable assumptions, it is possible to explain the variability in the observed range of braking indices of young pulsars
and to predict their evolutionary paths in the $P$-$\dot{P}$ diagram.

\subsection{Braking index and evolution timescale for realistic magnetic field evolution models.}

In Fig.~\ref{fig3} we plot the braking index evolution for three representative models,
Note that the horizontal axis in Fig.~\ref{fig1} represents the characteristic age $\tau_c$, while in Fig.~\ref{fig3} we show our results as a function of the real age of each model, so that a direct comparison is not possible.
We also used a linear scale in this plot to show the quasi-periodic oscillations during the long-term evolution more clearly.
For fields $\lesssim 10^{13}$ G and simple dipolar geometries (dashed line) the braking index at late times is
$n>3$, but its absolute value is low. In contrast, for strong dipolar fields, or for weak dipolar components but
with strong higher order multipoles, Hall-drift induced
oscillations appear sooner or later and,  in some situations, have large amplitudes that result in very high absolute values
of the braking index. It is particularly interesting to compare models A and C (dashed and solid line),  which have the same
initial dipolar component, representative of a typical pulsar (at $10^6$ yr the dipolar field is about $3$-$5\times 10^{12}$ G).
However, the presence of a strong octupolar component at birth results in a radically different braking index behaviour,
even if $P$ and $\dot{P}$ are similar.
The amplitude of the oscillations and whether the modes are damped or excited depends on particular details of the
small-scale structure of the magnetic field, which is unknown. The exact age at which these oscillatory modes are excited
is connected to the temperature of the star, and therefore to its cooling history and internal physics (neutrino emission
processes, superfluid gaps, etc.). For the standard cooling scenario this happens at $\approx 10^5$ yr. We can also
observe that as the star evolves and the magnetic diffusivity decreases, the frequency of the dominant
mode may vary (see solid line). For example, at early times in model C we observe
variability on shorter timescales (few kyr), while later we see longer oscillation periods ($10^5$ yr). After 1.5 Myr
higher frequency modes seem to be growing again.

\begin{figure}[th]
\centering
\includegraphics[width=0.45\textwidth]{fig4a.eps}
\includegraphics[width=0.45\textwidth]{fig4b.eps}
\caption{Phase residuals for model C around $t_0=10.93$ kyr (upper panel) and $t_0=100.5$ kyr (lower panel),
with the corresponding value of $n$ indicated in each case.
We show residuals after removing the best-fit solution including quadratic terms (solid) or including cubic terms (long dashes); dotted lines show the cubic residuals obtained  assuming a constant $B_0$.}
\label{fig4}
\end{figure}

\subsection{Timing residuals.}
In timing analysis of radio pulsars, an important piece of information is the study of the residuals, i.e, the phase differences between the observed signal and the best-fit model, including the frequency and frequency derivative. As shown in detail in \cite{Hobbs10}, there is a rich variety of shapes in the residuals of radiopulsars.
As we have shown, the strong Hall-induced interplay between different multipoles and the toroidal field produces a complex evolution of $B_0$. To compare with observational timing analysis, we
proceeded as follows: first, we choose a short interval of 40 yr in our simulations, centered on a fixed time denoted
by $t_0$. This interval is similar to the longest periods for
which phase-coherent timing analysis for radio pulsars can be performed.
Then, from our theoretical $B_0(t)$, we obtain $\nu(t)$
by integration of the classical spindown formula (Eq.~\ref{equ:spin-down-B(t)}) assuming an orthogonal rotator,
and the phase $\Phi(t)=\int_{t_0}^{t} \nu(t')dt'$. Finally, we fit our
synthetic time-dependent phase with a quadratic function
\begin{equation}
\Phi(t) = \Phi_0 + \nu_f (t-t_0) + \frac{1}{2}\dot{\nu}_f (t-t_0)^2 ~.
\end{equation}
The results of the fit, $\nu_f$ and $\dot{\nu}_f$, are of course in the range of values
of the {\it real} quantities in the time interval.

In Fig.~\ref{fig4} we show the phase residuals for model C (solid lines) at $t_0=10.93$ kyr (top) and $t_0=100.5$ kyr (bottom).
During the two analyzed intervals, the mean values of the dipolar magnetic field are $\bar{B_0}=8.1\times 10^{12}$ G and $\bar{B_0}=5.3\times 10^{12}$ G. The cubic shape
indicates that the residuals are dominated by the next term in the Taylor series (i.e., red noise). Including in the
fitting function a cubic term allows one to measure $\ddot{\nu}$ and the braking index ($n=0.39$ and $n=3.45$,
respectively). With dashes we show the fourth-order residuals after subtracting the third-order term in the fitting function.
For comparison, we also show with dotted lines the third-order residuals obtained assuming a constant value $B_0(t)=\bar{B_0}$ (that leads to $n=3$).

In the first time interval, the change of magnetic field, $\delta B_0=1.3\times 10^{10}$ G $=1.6\times 10^{-3} \bar{B_0}$,
is strong enough for the cubic residuals to be visibly different from the constant field case. The value of $|\ddot{\nu}|$
(and of $|n|$) is low, and the residuals of fourth order are significant. In the second time interval, $\delta B_0=1.4\times 10^8$ G $=2.6\times 10^{-5} \bar{B_0}$, and $n$ is close to 3. As a consequence, the deviation from the constant field spin-down behavior is slow and fourth-order residuals are orders of magnitude weaker (dashed line).

\begin{figure}[t]
\centering
\includegraphics[width=0.42\textwidth]{fig5.eps}
\caption{Comparison of the cubic residuals of the theoretical model (solid line) with the result after adding
a sinusoidal perturbation with frequency $0.1$ yr$^{-1}$ (dashed) or a random perturbation (dot-dashed).
For both perturbation forms we considered a maximum amplitude of $(\delta B_0/2)=6.5\times 10^9$ G.
}
\label{fig5}
\end{figure}

To investigate the effect of short-time irregularities, we repeated the process but
artificially added two additional sources of ``noise"  ($B_p$) to our theoretical values of $B_0(t)$,
with amplitude of the same order as the theoretical variations $\delta B_0$.
In Fig.~\ref{fig5} we show the residuals obtained by adding a sinusoidal perturbation,
$B_p=(\delta B_0/2)\sin(2\pi t/T)$, with $T=10$ yr (dashed lines), or a
random perturbation of maximum absolute value $\delta B_0/2$ (dash-dotted lines).
In our analysis, we sampled values every $\sim 10^{-2}$ yr, comparable with the typical integration time in observational data.
We checked that increasing the time interval strongly reduces
the contribution of the random noise at a fixed amplitude, as expected. Similarly, the periodic short-term noise
is reduced when the integration interval becomes much longer than the perturbation period.

Note that small perturbations of $B_0$ produce visible effects,
but the cubic residual still dominates. In constrast, in many pulsars a smooth cubic phase residual is clearly seen,
which means that the spin-down is very stable, short-term irregularities are completely negligible (or averaged out in
the long observation period), and the measure of a braking index (even a high value) dominated by the
secular evolution is robust. Glitches should appear as cusps in these plots (e.g. see data for PSR B0154+61
in \cite{Hobbs10}), but a more detailed investigation in this direction is beyond the scope of this paper.

\section{Conclusions.}

According to our current understanding of the magneto-thermal evolution of neutron stars, the time variation of the
magnetic field notably affects the braking index of pulsars and
contributes to the so-called red noise in the timing residuals. The qualitative picture seems to agree with the data:
  i) there is a short stage in which the dipolar field appears to be increasing, possibly reflecting the reemergence of magnetic field after initial accretion;
 ii) a stage dominated by Ohmic dissipation, while the star is still warm, and during which we always expect $n>3$;
iii) once the star has cooled down, after about $10^5$ yr, oscillatory Hall-drift modes in the crust are excited, and quasi-periodic oscillations of the braking index are consequently expected.

This is by no means the only contribution to the timing noise, nor necessarily the dominant one.
For example, changes in the effective moment of inertia as the superconducting region
of the core grows and unknown mechanisms for alignment of the rotation and magnetic field
\citep{TM1998,Young2010} would have similar long-term effects.
In addition to this secular evolution, short-term magnetospheric effects, such as oscillations of the inclination angle around the equilibrium position, or magnetospheric changes \citep{Spit2012L,Spit2012}, could be also responsible for the strong quasi-periodic features seen in phase residuals for several cases \citep{Lyne2010}.

Our results show that realistic magnetic field evolution models predict temporal
variations consistent with the observed trends. In particular, configurations similar to our model C could explain the
wide range in positive and negative values that are observed. This model has a moderate dipolar component
($10^{13}$ G at birth, which becomes a few times $10^{12}$ G at middle age), typical of pulsars, but a strong
octupolar component, large enough to activate Hall modes. The same generic behavior is expected for other models
with strong, smaller scale components, even if the dipole is weak.
However,  it must be mentioned that the models we studied
predict high values of the braking index ($10$-$100$) but not extremely high values ($10^3$-$10^4$), as reported
for pulsars with old characteristic ages. The absolute value of the braking index at late times, which can be
of either sign, is not simple to predict, since it depends on particular details of important components of the
magnetic field (higher order multipoles, toroidal component). Perhaps different field geometries
(e.g., without axial symmetry) and/or magnetospheric corrections may account for the discrepancy.

We emphasize that there is absolutely no reason to expect that the magnetic field remains constant during
a pulsar lifetime, and we must abandon the oversimplified models that assume constant (in time) pure dipolar fields.
On the other hand, since there is a strong interplay between the magnetic, thermal, and rotational evolution
of neutron stars, there are unexpected ways to obtain information. If future observations (or a deeper analysis of
archival data) increase the number
of reliable measures of braking indices, and we can firmly establish at which age the transition from the purely diffusive to the
oscillatory regime happens, there is coded information about the cooling history of neutron stars of
potential interest to place constraints on the interior physics.
In this respect, it is also important to mention that precise phase-connected
timing solutions can also be performed in other bands, as has been recently done for 17 $\gamma$-ray-selected pulsars
\citep{Ray2011}, using the Large Area Telescope (LAT) on the Fermi Gamma-ray Space Telescope. This and other
future high-energy missions can enlarge our data sample and improve the precision of our measurements.

\begin{acknowledgements}
This research was supported by the grants AYA 2010-21097-C03-02 and ACOMP/2012/135.
DV is supported by a fellowship from the \textit{Prometeo} program for research groups of excellence of the Generalitat Valenciana (Prometeo/2009/103). We thank C. Espinoza and J.A. Miralles for their valuable comments.
\end{acknowledgements}


\section{Introduction}

The recent use of three-dimensional Doppler tomography to investigate the Algol binaries U CrB \citepalias{agafonov09} and RS Vul \citepalias{richards10} has unveiled a significant amount of detail regarding the accretion process in these systems. As well as confirming the presence of a number of features seen in two-dimensional tomographic data, the new data have shown an unexpectedly large amount of gas located outside of the orbital plane in both systems. U CrB also displays a transient accretion disc inclined with respect to the orbital plane, which is not expected if mass transfer occurs as a result of gravitational and centrifugal forces alone. The companions in short-period Algols are rapidly rotating late-type stars that are likely to be magnetically active.  Magnetic behavior occuring at or near the first Lagrange point (L1) may offer an explanation for the out-of-plane material.

The notion that large-scale magnetic effects may drastically influence the morphology and time-evolution of the accretion processes is evidenced by the observations of \citet{richards10}, who noted the similarity in the out-of-plane velocities located at L1 and magnetically active regions on the donor star in RS Vul. Moreover, they found evidence of a loop prominence on the cool star close to the L1 point in the 3D tomogram of RS Vul, which suggests that the magnetic field may have deflected the gas stream relative to the central plane. Such a deflection would cause gas to be introduced into the system with additional (non-z) angular momentum that could then propagate toward the primary and persist as a tilted disc. This deflection could be attributed to a variety of mechanisms. Starspots on the donor in the prototype Algol system $\beta$ Persei were observed by \citet{olson81}, and coronal ejections in eight other Algol systems, including RS Vul, have been detected by \citet{white83}. More recently, a large-scale persistent coronal loop in Algol was fully resolved by \citet{peterson10}, suggesting the possibility of an extended magnetic field structure between the two stars. While the likelihood of magnetic activity in U CrB and RS Vul is high, the lack of information presently available regarding the exact magnetic structure of these systems makes it difficult to determine the degree to which magnetic effects are responsible for the transport of gas out of the orbital plane.

Further complicating the mass transfer process in Algols is the fact that their short orbital period ($< 5-6$ days) combined with the close orbit and relatively large radius of the primary results in the accretion stream directly impacting or grazing the surface of the accreting star. Direct-impact systems display a permanent or transient Keplerian disc and sometimes a partial disc or localized region characterized by an asymmetric mass distribution and sub-Keplerian velocities \citep{richards92solo, richards99}. In these systems, the localized region is associated with gas that has circled the accreting star and has slowed down because of the interaction with the incoming gas stream \citep{richards92solo}. Gas is expected to leave the orbital plane due to interactions at the site where the stream impacts the surface of the primary, the stream-disc interaction, and the localized region \citep{richards04}.

Both U CrB and RS Vul are short-period binaries, and both display a number of H$\alpha$ emission features characteristic of direct-impact systems. U CrB has been classified as an alternating system, since it has been observed in both a streamlike accretion state and a disclike state \citep{richards99}. The stream state is characterized by strong emission along the predicted ballistic trajectory of the accretion stream against a lack of circumprimary emission. In the disc mode, emission from the accretion disc is the most prominent feature. The mechanism responsible for the transition between these two states has yet to be explained, although hydrodynamic simulations of $\beta$ Per \citepalias{blondin95} showed that the morphology of the accretion flow depended strongly on the density of the inflowing gas. The hydrodynamic simulations of \citet{richards98} also measured a large relative increase in the H$\alpha$ emissivity of the stream relative to the disc when stream density was varied. The similarity of the stellar components and orbital parameters of RS Vul to those of U CrB (RS Vul = B5V + G1 III-IV, P=4.48 days; U CrB = B6V + G0 III-IV, P=3.45 days) makes it likely that a disc mode is possible in RS Vul, and it has been classified as an alternating Algol by \citet{richards99} despite having not yet been observed in a disc mode.

Three-dimensional tomographic imaging of U CrB has allowed for a more thorough analysis of the accretion features by enabling the creation of 3D velocity space plots of the H$\alpha$ emission sources. Mapping this data to Cartesian coordinates is not possible due to the inability to resolve the exact Cartesian structure of the system, but the theoretical analysis of key features of the system provides reference points from which the gas dynamics of the system can be inferred. Features commonly used as landmarks include the motion of the stellar components, the ballistic trajectory of the accretion stream, and the locus of the accretion disc as defined by the Keplerian velocities at the surface of the primary and its outermost Roche surface.

The first 3D tomographic analysis of U CrB was performed by \citetalias{agafonov06}, and showed the system in a disc mode. An additional series of H$\alpha$ profiles were later analyzed in three dimensions and displayed a stream mode \citep{agafonov09}. Notable features of these two modes include: (1) A gas stream transferring mass from the donor to the primary star with velocities as high as 450 km s$^{-1}$ in the orbital plane and negligible z-velocities. (2) Circumprimary emission associated with an accretion disc that displays both positive and negative z-velocities ($v_z$ = $\pm 420$ km s$^{-1}$). The different signs of the z-velocities on opposite sides of the primary star suggest that the disc may be inclined with respect to the orbital plane. (3) A stream-disc interaction region, or ``hotspot'' ($v_z$ = +180 km s$^{-1}$ to +540 km s$^{-1}$), as well as a star-stream impact site ($v_z$ = -300 km s$^{-1}$ to +60 km s$^{-1}$). (4) A localized region ($v_z$ = +200 km s$^{-1}$ to +540 km s$^{-1}$). In contrast to the two accretion modes displayed by U CrB, RS Vul has only been observed in a stream state with a much stronger out-of-plane flow. A three-dimensional tomographic analysis of RS Vul was performed by \citet{richards10}, and displays many of the same features as U CrB. Major differences from U CrB include: a truncated gas stream ($v_z = -85$ km s$^{-1}$) with no indication of a star-stream hotspot, asymmetric circumprimary emission with no well-defined disc, a localized region diametrically opposite the circumprimary emission with opposite z-velocities (suggesting that the gas orbiting the accreting star is flowing along a tilted path with $v_z = +15$ to $-85$ km s$^{-1}$), and emission from a magnetically active region on the donor that suggests the presence of a loop prominence ($v_z = -150$ to $150$ km s$^{-1}$).

The fact that the H$\alpha$ emission at L1 possesses z-velocities equal to those of magnetically active regions on the donor strongly suggests that the stream is magnetically deflected at or near L1. The relatively high velocities associated with the stream also suggest that the stream velocity may be boosted to many times the speed that would occur via Roche lobe overflow alone. Our goal was to determine whether a deflected stream is responsible for introducing a sufficient amount of tilt to the accretion disc, and to what extent the deflection influences the presence or absence of other out-of-plane accretion features. To investigate this we performed hydrodynamic simulations of mass transfer in RS Vul in which magnetic effects were mimicked by manually deflecting and boosting the stream at L1. The resulting steady-state solutions were analyzed to determine the effect of the deflection on the morphology of the circumstellar material, the presence of out-of-plane flow, and the presence or absence of the accretion features noted in the tomographic data. The hydrodynamic model we used is described in section 2, results from runs with undeflected and deflected streams are given in section 3, an analysis of the tilt of the disc with comparisons to analytical estimates is given in section 4, and our results summarized in section 5.

\section{Computational Methods}

\subsection{Grid parameters}

Our simulations were performed with the three-dimensional hydrodynamics code VH-1, which was developed by the Numerical Astrophysics Group at the Virginia Institute for Theoretical Astrophysics. VH-1 solves the Euler equations for the flow of an inviscid fluid characterized by an adiabatic index $\gamma$. To solve the system, VH-1 utilizes operator splitting to reduce the grid into a number of 1-D segments. The fluid equations are solved along each segment using the Piecewise Parabolic Method, an explicit high-order Godunov method employed on an Eulerian (fixed) grid \citep{colella84}. The PPM can be used for both 2D and 3D geometries without necessitating significant changes to the algorithm other than the extension from 2D polar to 3D spherical coordinates.

We adapted VH-1 to model our system in spherical coordinates on a corotating grid located in the reference frame of the binary system and centered upon the primary star.  The inner boundary of the grid coincides with the surface of the primary at a radius of $R_1$, and the outer boundary extends to twice the distance from the centre of the primary to L1 ($\sim 6.2 R_1$). Since the system does not vary widely in scale between the two stars, we space our zones equally in the radial direction. The parameters of U CrB and RS Vul are given in Table 1, and were taken from \citet{agafonov09} and \citet{richards10}. We note that while we use the parameters corresponding to RS Vul in our simulations, both the orbital parameters and stellar components of U CrB are similar to those of RS Vul. This enables us to interpret our results in the context of either system.

To avoid grid convergence and the coordinate singularities intrinsic to spherical coordinates, we employ a yin-yang overset grid \citep{kageyama04}. The yin-yang grid consists of two separate spherical polar grids oriented 90$^\circ$ with respect to each other, so that the polar region of one grid aligns with the equatorial region of the other. Each grid extends $\pm \pi/4$ azimuthally, and $\pm 3\pi/4$ in the polar direction (Figure \ref{fig:yinyang}). This grid configuration allows for reasonable limits on the Courant number and ensures that our runtimes do not become excessive. As a result, we are able to extend our outer boundary to a sufficiently large radius without imposing limits on the the azimuthal extent of the grid. We define our coordinates such that the positive x-axis of the yin grid lies along the binary axis and points toward the secondary. The positive y-axis points in the direction of motion of the secondary, and the z-axis is perpendicular to the orbital plane. Each piece of the grid has dimensions of $288 \times 144 \times 432$, which results in an angular resolution of $\Delta\phi \sim 0.625^\circ$ and $\Delta r/R_1$ = 0.018 at the inner boundary.

\begin{figure}
\includegraphics[width=84mm]{yinyang}
\caption{Yin-yang grid geometry. The figure on the left shows the range of one grid, and the figure on the right shows the overlap between the two grids.}
\label{fig:yinyang}
\end{figure}

Since the inner boundary corresponds to the surface of the accreting star, we must choose a boundary condition that allows gas to accrete while still adequately modeling the effects of the stellar atmosphere. The fact that the stream penetrates only a small portion of the stellar atmosphere \citep[][$\sim 10^{-6}R_1$ for Algol, ]{richards92solo} allows us to approximate the surface of the primary as a solid surface. Interactions between the atmosphere and the stream could have a non-negligible effect on the structure of the circumstellar gas (especially in the case of a rotating and precessing primary), but for the purposes of this paper we neglect the structure of the atmosphere. Instead, we employ a partially reflecting boundary condition that sets the radial velocity of the gas at the innermost radial zone to zero every timestep. This approximates the gas encountering a shock at the surface of the star and accreting shortly afterward. At the outer boundary, gas should be allowed to freely exit the grid, so we implement a modified zero-gradient boundary condition with the added requirement that gas may not be allowed to flow inward once it has passed off the grid.

\begin{table*}
\begin{minipage}{126mm}
\caption{Parameters for U CrB and RS Vul}
\label{symbols}
\begin{tabular}{lll}
\hline
Parameter & U CrB & RS Vul \\
\hline
Mass of Primary, $M_1$ & $4.80 M_{\odot}$ & $6.59 M_{\odot}$ \\
Mass of Secondary, $M_2$ & $1.40 M_{\odot} $ & $1.76 M_{\odot}$ \\
Effective Temperature of Primary, $T_1$ & $14000$ K & $16500$ K \\
Effective Temperature of Secondary, $T_2$ & $5250$ K & $6000$ K \\
Radius of Primary, $R_1$ & $3.0 R_{\odot}$ & $4.7 R_{\odot}$ \\
Radius of Secondary, $R_2$ & $4.6 R_{\odot} $ & $5.8 R_{\odot}$ \\
Binary Separation, $a$ & $1.23 \times 10^{12}$ cm & $1.61 \times 10^{12}$ cm \\
Orbital Period, $P_{orb}$ & $3.452$ days & $4.478$ days \\
\hline
\end{tabular}
\end{minipage}
\end{table*}

\subsection{Radiative cooling and stream density}

The 2D simulations of Algol by \citet{blondin95} found that the morphology of the circumprimary flow depended strongly on the initial stream density due to the effects of radiative cooling. The stream density's influence on the flow is seen primarily in the stream-disc interaction. A denser stream will cool more rapidly at the shock, negating the effects of heating within the shock. This results in a flow that behaves similarly to a compressible isothermal gas. A lower-density stream will not cool to the same extent, and will behave adiabatically.

To investigate the relative effect of cooling in U CrB and RS Vul, we compare the cooling time at the shock to a characteristic dynamical time. The dynamical time can be found using the gas speed at the outer edge of the disc where the stream-disc shock is located. Assuming a Keplerian velocity at the outer edge of the disc ($\approx 350$ km s$^-1$) allows us to calculate a dynamical time of:

\begin{equation}
t_d = \frac{GM_1}{v_{Kep}^3}
\end{equation}

In both U CrB and RS Vul, the dynamic timescale is calculated to be $0.04 P_{orb}$.

If a cooling function is known, we can calculate the characteristic cooling time at the shock to be

\begin{equation}
t_c = \frac{kT}{n_s\Lambda(T)}
\end{equation}

Here, $n_s$ is the electron number density in the post-shock flow, and $\Lambda$ is the cooling function evaluated at the temperature of the shocked gas, $T$. We calculate $\Lambda(T)$ with the cooling function for a low-density plasma as given by \citet{cox71}. From the Rankine-Hugoniot jump conditions, we can calculate the post-shock temperature to be $T=3\tilde{m}v^2_{sh}/16k$, where $\tilde{m}$ is the average mass per particle. We assume a gas comprised largely of ionized Hydrogen, in which case $\tilde{m} \sim m_H/2$. We assume the velocity of the shocked gas, $v_{sh}$ to be on the order of the Keplerian velocity at the shock site. For both U CrB and RS Vul, we calculate $T_{sh}$ to be approximately $1 \times 10^6$ K, which results in a cooling coefficient of $\Lambda = 8 \times 10^{-23}$ ergs cm$^3$ s$^{-1}$. An analytical estimate of the density at the shock site is difficult to obtain, so we note that our simulations show densities in the shock region on the order of $\rho = 0.1 \rho_{L1}$. Converting this density to physical units requires a measurement of the mass transfer rate. Since data on the exact mass transfer rates for U CrB and RS Vul are not available, we use the lower limit of the transfer rate in Algol, $\dot{M} \approx 10^{-11} M_{\odot}$ yr$^{-1}$ \citep{richards92solo}.

Comparing the characteristic dynamical time to the cooling time, we find that $t_d / t_c = 2.4$. This is a conservative estimate, as we have used a small value of the mass transfer rate to calculate the cooling time. Any increase in the mass transfer rate will serve to decrease the cooling time further, since $t_c$ varies as $1/\dot{M}$. In general, the mass transfer rates in short-period Algols vary widely between systems, and are not likely to remain constant in time. For Algol, $\dot{M}$ could range from $10^{-11}$ to $10^{-9} M_{\odot}$ yr$^{-1}$ \citep{richards92solo}. With these factors taken into account, it is likely that cooling at the shock is sufficiently rapid to counter the heating that occurs in the stream-disc interaction. This rapid cooling ensures that the shock-heated gas will quickly fall to some minimum temperature below $10^4$ K, at which point the recombination of H will diminish the cooling rate. The cooled gas will be heated by the strong photoionizing UV flux from the primary, which is sufficient to maintain the temperature at a value close to $10^4$ K. Taking this information into account, we assume a constant temperature of $11,000$K  in our simulations. This value falls midway between the effective temperatures of the primary and secondary in both U CrB and RS Vul. This allows us to expedite our calculations by employing an isothermal equation of state ($\gamma = 1$).

\subsection{Initial conditions at L1}

\citet[hereafter LS]{lubow75} conducted an extensive analysis of the gas dynamics at the L1 point, and we use their prescription to determine the initial conditions of the accretion stream. LS parametrized the flow in terms of the ratio of isothermal sound speed to the system's orbital speed: $\varepsilon = c_s / \omega a = (kT/m)^{1/2}$, for which we obtain $\varepsilon=0.035$. This gives a local sound speed of $9.5$ km s$^{-1}$ at L1. LS calculate that the stream will be accelerated to supersonic speeds as it travels through L1, so we assume an unboosted speed of Mach 1. The angle the stream makes with the binary axis can also be calculated using $\varepsilon$, and for RS Vul we obtain $\theta = 20.4^{\circ}$. Assuming an ideal gas, the pressure at L1 should be initialized to a value on the order of $\rho \varepsilon^2$.

For this work we require control over the direction and magnitude of the gas velocity at L1 in order to manually deflect the stream. To achieve this, we calculate the position of L1 on the grid and hold the stream parameters in the corresponding zones fixed at the values given above. This eliminates the need for the donor star to be modeled hydrodynamically, and we remove it from the grid while retaining its term in the gravitational potential. A small amount of flow into the region corresponding to the secondary occurs, but it remains confined to the volume that the secondary would normally occupy and does not affect the behavior of the accretion flow. In addition to giving us more control over the stream at L1, this strategy has the benefit of allowing for a more economic grid size, since the full secondary need not be present on the grid.

LS calculate the width of the stream at L1 to be on the order of $\varepsilon a$, which requires an angular resolution of at least $\Delta\phi = 0.022$. The angular resolution of our grid is small enough to accomodate this with two zones; each has an angular extent of $\Delta\phi = 0.011$. Symmetry requires the gas stream at L1 to extend equal distances above and below the orbital plane as well as an equal distance to either side of the binary axis. This necessitates the use of four zones to represent L1. In reality the density profile of the gas at L1 is Gaussian, but we lack the resolution to resolve this. Instead we hold the density of the gas fixed to a constant at our L1 zones. Ultimately, our system exhibits behavior similar to that of \citet{blondin95} in that the morphology of the accretion flow is relatively insensitive to the geometry of the zones at L1.

Our simulations do not assume any gas to be initially present around the primary, and the system requires many orbital periods to allow the full formation of a disc. A steady-state solution in which a disc is present and has stopped gaining mass (indicated by $\dot{M}_{acc} = \dot{M}_{L1}$) requires $40 - 50$ P$_{orb}$ according to preliminary low-resolution runs. To expedite the disc formation, we multiply the density in the disc by a constant factor after two orbital periods have elapsed. Since the mass accretion rate is proportional to the density at the inner boundary, this results in the mass accretion rate being boosted to a value closer to that at which it stabilizes.

\section{Results}

\subsection{Undeflected stream}

To ensure that our model behaved appropriately in the absence of a deflected stream, we performed a fiducial run with the stream left unmodified ($\phi=0^{\circ}$, $v_{stream}=1M$). Figure \ref{fig:notilt_topside} shows the stream following the ballistic trajectory as it falls toward the primary and narrowing as it approaches the surface. The stream impacts the surface after approximately $0.1P_{orb}$. We allow gas to accrete until $2 P_{orb}$, during which a thin accretion disc begins to form. At $2 P_{orb}$, we boost the disc density. The boost instigates some radial density oscillations due to increased pressure along the density discontinuity that occurs at the radius of the boosted region. These oscillations propagate with the flow around the disc, but fade within the next few orbital periods. Following the density boost, the radius of the disc continues to increase for several orbital periods, and a steady state is achieved at roughly $15 P_{orb}$.

\begin{figure*}
\begin{minipage}{126mm}
\includegraphics[width=126mm]{notilt_topside_v3b}
\caption{Top and side views of the primary star and and disc for an undeflected stream after $18 P_{orb}$. The density scale is normalized to 1 at L1. The top view depicts the density on a slice through the orbital plane. Density isosurfaces in the side view are taken at $\rho = 1 \times 10^{-1}, 5 \times 10^{-2}, 1 \times 10^{-2}, 5 \times 10^{-3},$ and $1 \times 10^{-3}$ with the transparency increased by a factor of 0.5 for each isosurface.}
\label{fig:notilt_topside}
\end{minipage}
\end{figure*}

When the system has stabilized, we measure the disc to have an outer radius of $0.34a$ (where $a$ is the binary separation as measured from the center of the accreting star). Examining the density profile averaged over $\phi$ (Figure \ref{fig:density_v_radius}) shows that the density of the innermost portion of the disc ($\sim 0.22a$) is of the same order of magnitude as the density at L1 ($\sim 0.63a$). The density decreases sharply to just above $\sim 10^{-4} \rho_{L1}$ at the outer edge of the disc. The disc is symmetric about the orbital plane to approximately $\pm 0.03a$ above and below the plane. Beyond this region the density decreases rapidly, and is measured to be less than $10^{-4}\rho_{L1}$.

The velocities within the disc are oriented largely within the plane. As the gas accelerates from L1, its speed increases from $\sim 10$ km s$^{-1}$ to a final velocity of $\sim 240$ km s$^{-1}$ before encountering the stream-disc shock. (All velocities are measured with respect to a corotating frame.) After crossing the shock and entering the higher-density regions of the disc, the gas assumes velocities which fall within the range of $240-440$ km s$^{-1}$, with small velocities corresponding to the outer edge of the disc, and higher velocities located near the surface of the accretor. The velocities within the disc remain relatively constant in direction and are sub-Keplerian throughout ($60- 80$ percent at the outer and inner edges of the disc, respectively). These outer velocities are less uniform near the localized region and stream-disc interaction site, and reach values up to $300$ km s$^{-1}$ in the vicinity of the stream-disc shock.

\begin{figure}
\includegraphics[width=84mm]{density_v_radius}
\caption{Density vs. radius, averaged over $\phi$ within the orbital plane at $18 P_{orb}$. Densities are normalized to 1 at L1. The solid line depicts the density for the undeflected stream, while the dashed line corresponds to a stream tilted to $45^{\circ}$ and boosted to Mach 30. L1 is located at $0.63a$, and the surface of the primary at $0.20a$.}
\label{fig:density_v_radius}
\end{figure}

To quantify the angle with which the disc is inclined, we integrate the angular momentum present in the disc and measure its angle with respect to the z-axis. We limit our integration to the region occupied by the high-density portions of the disc, using Figure \ref{fig:density_v_radius} as a guideline. The result is summed angular momentum values for each value of $(r, \phi)$ in the orbital plane. Each value of $r$ defines a ring within the disc for which we are able to calculate the tilt. The sign of the tilt is set equal to the sign of the x-component of the net angular momentum. This has the effect of assigning positive tilt angles to rings that are tilted toward the donor, and negative angles for rings that are tilted away. For the undeflected stream, the disc is largely symmetric about the orbital plane and remains essentially untilted over time. Examining the tilt as a function of radius at $18 P_{orb}$ (Figure \ref{fig:tilt_v_radius}) shows that the high-density portions of the disc from the inner boundary to $\sim 0.3a$ only possess tilt on the order of $0.1^{\circ}$. This is a fraction of the angular resolution of by one grid zone. The tilt  varies gradually as the radius increases, and the sign changes three times along the radius, indicating that some warping occurs about the orbital plane. The tilt becomes slightly larger in the low-density portions of the disc and reaches a maximum value of $0.5^{\circ}$. This tilt is smaller than the amount able to be resolved with one zone. The presence of this tilt is likely a result of a comparatively large amount of non-z angular momentum created in the stream-disc interaction through the splashing of low-density gas. Turbulence in the gas flowing along the outer low-density edge of the disc may also contribute to this measurement.

\begin{figure}
\includegraphics[width=84mm]{tilt_v_radius}
\caption{Tilt of the disc vs. radius, calculated at $18 P_{orb}$ by summing the angular momentum within the disc over $\phi$ for each value of $r$. The solid line depicts tilt for the undeflected stream, while the dashed line corresponds to a stream tilted to $45^{\circ}$ and boosted to Mach 30.}
\label{fig:tilt_v_radius}
\end{figure}

The disc's orientation with respect to the orbital plane can be found by examining a plot of the x and y components of the angular momentum. In cartesian space, the projection of the angular momentum vector into the orbital plane points toward the direction that the disc is azimuthally inclined. (This is analogous to treating the disc as a rigid body with a rod through its axis of rotation. The projection of the angular momentum vector into the orbital plane is the direction you would push the rod to tilt the disc.) Figure \ref{fig:orientation} shows a $J_y$ vs. $J_x$ plot for both runs, with $J_x$ and $J_y$ components shown as a fraction of the total angular momentum. The data points for the undeflected stream lie clustered around the origin, with magnitudes of less than $1$ percent of the total angular momentum. This supports the qualitative observation that the disc lies almost entirely in the orbital plane throughout the length of the simulation. The disc appears to precess at a period of $P_{orb}$, which is a by-product of measurements of the angular momentum taken from within a co-rotating reference frame.

\begin{figure}
\includegraphics[width=84mm]{orientation}
\caption{Components of the angular momentum of the disc for the undeflected stream (circles) and the deflected and boosted stream (crosses). Data points were taken at intervals of $\sim 0.2 P_{orb}$. The undeflected stream results in an disc that remains largely in the orbital plane, with minor variations over time that remain symmetric about the origin. The deflected and boosted stream produces a disc biased toward $+J_x$ due to the stream impacting the far side of the disc from below and tilting it toward the donor star.}
\label{fig:orientation}
\end{figure}

\subsection{Deflected stream}

To simulate a system with a deflected and a boosted velocity, we modified the stream parameters at the zones corresponding to L1 by directing the stream velocity $45^{\circ}$ below the orbital plane and boosting it to thirty times the local sound speed; this is comparable to gas flows within magnetic loops. These values were chosen to maximize the likelihood that a tilted disc would form. In the early stages of evolution (Figure \ref{fig:tilt_topside}), the stream falls along the predicted ballistic trajectory, which extends below the orbital plane to just below the bottom of the accreting star. As it nears the primary, it is pulled back toward the orbital plane and wraps around the accretor. The boosted velocity is large enough to cause the stream to miss the surface of the star, resulting in the absence of a star-stream collision and a gas stream that extends above the orbital plane to nearly $2 a$. The stream then circles the accretor and interferes with the gas stream from the secondary as an accretion disc begins to form near the orbital plane. The flow during the first few orbital periods displays a large degree of out-of-plane flow as the system approaches a steady state. At $2 P_{orb}$ the density within the disc is boosted, inducing some minor radial oscillations as in the undeflected case. As the disc gains mass, the stream-disc interaction continues to drive gas upward above the orbital plane. The system stabilizes over the next $15 P_{orb}$, after which the accretion disc remains present along with a large plume of gas deflected above the orbital plane by the stream-disc interaction.

\begin{figure*}
\begin{minipage}{126mm}
\includegraphics[width=126mm]{tilt_topside_v3b}
\caption{Top and side views of the primary star and disc for a stream boosted to Mach 30 and deflected $45^{\circ}$ below the orbital plane. The density scale is normalized to 1 at L1. Density isosurfaces are taken at $\rho = 1 \times 10^{-1}, 5 \times 10^{-2}, 1 \times 10^{-2}, 5 \times 10^{-3}, 1 \times 10^{-3}$ with the transparency increased by a factor of 0.5 for each isosurface. Due to asymmetry about the orbital plane, density isosurfaces are depicted in the top view, as opposed to a slice along the plane.}
\label{fig:tilt_topside}
\end{minipage}
\end{figure*}

Investigating the morphology of the accretion disc shows some radial enlargement and warping due to the deflected stream. Averaging the density in the central plane over $\phi$ (Figure \ref{fig:density_v_radius}) shows that the density falls off more slowly than in the undeflected case, reaching $10^{-3}\rho_{L1}$ at $0.5a$. The inner density of the disc remains on the same order of magnitude as in the undeflected case, and is slightly larger in radius. The top view of Figure \ref{fig:tilt_topside} at $18 P_{orb}$ shows that the disc is also asymmetric, with a slightly enlarged portion in the upper right portion of the disc. This large density region originates as the disc forms and propagates around the disc with time. It does not appear to originate due to the artificial density boost that occurs at $2 P_{orb}$. The vertical extent of the disc is similar to that of the undeflected case, neglecting slight enlargement above and below the orbital plane due to warping that occurs in the outer edges of the disc. The majority of the disc still lies in the orbital plane, but with lower density gas ($\sim 0.1 \rho_{L1}$) deflected above the plane due to the stream-disc interaction. This gas is pulled back into the plane, where further dissipation of angular momentum due to interactions with the disc causes a portion of the material to join the disc and eventually accrete.

The deflected stream begins at a speed of $285$ km s$^{-1}$ (Mach 30). As the circumprimary material begins to build up into a disc, gas is blown upward in the collision with z-velocities of up to $210$ km s$^{-1}$. Within the disc we observe a circular orbit with velocities that remain sub-Keplerian in the same range as the undeflected case ($240-440$ km s$^{-1}$ for the high-density portion, $150-250$ km s$^{-1}$ for the low-density portion). We also measure z-velocities on the order of $\pm $ 10 km s$^{-1}$ due to the warping and tilting in the disc.

Examining the tilt within the disc as a function of radius after $18 P_{orb}$ has elapsed (Figure \ref{fig:tilt_v_radius}) shows that the warping within the disc is more pronounced than in the undeflected stream. The inner portion of the disc is tilted at $\sim 1.5^{\circ}$ near the the surface of the primary, and increases to $\sim 3^{\circ}$ at $0.35a$. Warping is indicated by the variations in tilt angle with increasing radius. The tilt continues to increase at radii greater than $0.35a$, mostly due to the high z-velocities in the material blown upward in the stream-disc interaction. The tilt angle calculated for this low-density gas reaches values as high as $50^{\circ}$ at the outer edge of the disc. All the tilt measurements for this timestep result in positive values, indicating that the disc is tilted consistently in the $+x$ direction (toward the secondary).

The orientation of the disc changes noticeably over time when we deflect and boost the stream. The plot of $J_y$ vs. $J_x$ (Figure \ref{fig:orientation}) shows a much wider spread of values than that of the undeflected stream. The fraction of angular momentum that is distributed in the $x$ and $y$ components remains between 2 and 3 percent at all times, and grows to 5 percent when the disc is oriented toward the donor star. When this occurs, the effects of the stream's impact on the disc are enhanced due to the upward force of the stream on the already upwardly-moving gas in the disc. As in the undeflected case, the apparent precession is due to measurements from within the corotating reference frame.

\subsection{Emissivity}

Additional information regarding the out-of-plane velocities and accretion features are more easily visible if we mimic the tomographic results with a pseudo-emissivity plot in velocity space. This was first done for Algol systems by \citet{richards98}. In our simulations, emissivity scales as $\rho^2,$ and lacks any dependence on the temperature or ionization state of the gas. We slice our velocity space plots into slabs along $v_z$, sum the emissivity along $v_z$ within each slab, and scale the emissivity to the maximum value observed in the slab. This maximum value is also listed in terms of the maximum emissivity on the entire grid. The velocities measured in our hydrodynamic simulations take place in a corotating frame, and we transform these into an inertial frame for the sake of comparison with tomographic data. Since we use the inner boundary to represent the primary and we have removed the secondary from the grid, emission features associated with stellar activity occuring on the surface of either star (e.g. circumprimary emission, chromospheric emission from the secondary) are absent from the plot. As a result, the plots appear more sparse than the original tomographic data presented in \citet{agafonov09} and \citet{richards10}. We also note that the strong density dependence causes the emissivity in the central plane to appear many orders of magnitude greater than that of other slices.

The emissivity plots for the undeflected stream (Figure \ref{fig:halpha_nt}) clearly show the presence of the disc in the central slice ($v_z = -30$ to $30$ km s$^{-1}$). Since the inner portion of the disc circles the primary with a higher velocity than the outer portion, the disc is turned inside-out on the velocity space plot. The bulk of the disc possesses velocities that fall within the Keplerian predictions (400-515 km s$^{-1}$), although we also measure sub-Keplerian velocities as small as 370 km s$^{-1}$. The accretion stream is also visible in the central slice, beginning near $v_x = 0$ km s$^{-1}$ and extending further in the -x direction as it gains speed and falls along the binary axis. The slightly stronger emissivity region located along the inside region of the``ring'' at roughly $(v_x, v_y)$ = ($-250,250$) km s$^{-1}$ corresponds to the stream-disc interaction. Emission from this feature extends up to $v_z=30$ km s$^{-1}$. We also note that the emissivity is slightly enhanced near ($\pm 350,-350$) km s$^{-1}$ due to the overlap between the yin and yang grids.

Some emission is present in the slab corresponding to $v_z$ = 30 to 90 km s$^{-1}$, and is located primarily within the locus of the accretion disc. The intensities of this emission are seven orders of magnitude below the maximum intensity on the grid, and distinct emission features are not immediately apparent in this slice. For gas wtih downward velocities of $v_z$ = $-30$ to $-90$ km s$^{-1}$, some emission is present in the locus of the disc again, though at intensities ten orders of magnitude below the maximum grid intensity. For velocities greater than $v_z$ = $\pm$90 km s$^{-1}$, the emission lacks any distinct features and occurs at intensities on the order of $10^{-13}-10^{-14}$ of the maximum intensity.

\begin{figure*}
\begin{minipage}{126mm}
\includegraphics[width=126mm]{ev_rsvul1100_NT}
\caption{Emissivity for RS Vul with an undeflected stream after $18 P_{orb}$. The emission intensity on each plot is normalized to the maximum within the individual slice. The value of each maximum is also listed in terms of the maximum intensity on the entire grid.}
\label{fig:halpha_nt}
\end{minipage}
\end{figure*}

The emissivity plots of the tilted disc (Figure \ref{fig:halpha_tilt}) show much higher emissivity measurements in slices above and below $v_z=0$. The disc is again the dominant feature in the central slice, and is larger in radius than in the undeflected case, as expected from the density profiles. This results in an emissivity region that extends inward to a slightly greater extent. The majority of the velocities lie within the Keplerian predictions ($375-515$ km s$^{-1}$), with some gas possessing velocities as low as 350 km s$^{-1}$. Since the stream has been deflected below the orbital plane, it appears in the upper left quadrant of the $-v_z$ slices. The stream-disc interaction is visible in the slices corresponding to positive z-velocities, and appears most prominently near $(v_x, v_y)$ = ($-450,150$) km s$^{-1}$ in the slice of $v_z = 30$ to $90$ km s$^{-1}$. The high $v_z$ material on the far side of the disc is also apparent in the emissivity plots. As it rises and wraps around the secondary (with $+v_z$), it can be seen in the lower left quadrant. As it falls back towards the orbital plane ($-v_z$), it shifts toward the lower right quadrant, indicating motion back toward L1. The gas collides with the orbital plane again, and a portion is driven below the orbital plane once again. This is visible in the $+v_z$ slices as a feature in the upper right quadrant.

As a whole, the emissivity for the deflected and boosted stream reveals values above and below $v_z=0$ with intensities much closer to the values measured in the central slice. This indicates a much stronger out-of-plane flow produced by the deflected and boosted stream. It should also be noted that the maximum intensity on the grid is nearly five times larger for the deflected and boosted stream than for the undeflected stream. This is to be expected, since the densities within the disc are slightly larger when the stream was deflected, and any changes in density will appear in the emissivity as $\rho^2$. Examining the emissivity in the early stages of our run ($2 P_{orb}$, just before the density is boosted) shows the intensity more uniformly distributed throughout the $v_z$ slices. No notable differences in the structure of the emissivity features are seen.

\begin{figure*}
\begin{minipage}{126mm}
\includegraphics[width=126mm]{ev_rsvul1100_T}
\caption{Emissivity for RS Vul with an accretion stream deflected below the orbital plane by $45^{\circ}$ and boosted to Mach 30 at $18 P_{orb}$. As in the previous figure, the emission intensity is normalized to the maximum value within each individual slice, and the value of each maximum is listed in terms of the maximum intensity on the entire grid.}
\label{fig:halpha_tilt}
\end{minipage}
\end{figure*}

\section{Analysis}

The main features visible on our emissivity plots for both the undeflected and deflected streams are the disc emission, accretion stream, and stream-disc interaction. In the undeflected case we also measure non-localized ``clouds'' with z-velocities up to $\pm 210$ km s$^{-1}$, but with magnitudes smaller than the emission features located in the central plane by a minimum of eight orders of magnitude. For the deflected stream, the main sources of emissivity are the disc and the out-of-plane gas located on the far side of the secondary (extending from $v_z=+210$ km s$^{-1}$ to $v_z=-150$ km s$^{-1}$). The disc structure is similar to that of the undeflected case, but displays some radial enlargement and slight warping. The out-of-plane gas appears with much stronger intensities relative to the gas in the orbital plane, and measures $10^{-3}-10^{-7}$, as opposed to $10^{-8}-10^{-14}$.

The early stages of our deflected stream simulation are comparable to the tomographic data of RS Vul \citep{richards10}. The emission seen in the tomographic data of RS Vul is not confined to the Keplerian locus and is found primarily along the stream, near the mass donor and accretor, and in several other localized high-velocity regions. In our simulations, prior to the formation of the accretion disc (i.e. before the density boost occurs at $\sim 2 P_{orb}$) we observe a similar structure. The stream is visible in slices with negative velocities as low as $v_z=-210$ km s$^{-1}$, and the high-velocity features corresponding to the stream material wrapping around the primary are visible from $v_z = -150$ km s$^{-1}$ to $v_z = +210$ km s$^{-1}$. As in the tomographic data, the star-stream interaction is not visible, since the deflected stream misses the star entirely. Emission from the Keplerian locus is non-uniform and occurs at intensities lower than those corresponding to the stream by an order of magnitude.

Despite similarities with the morphology of our simulations, the tomographic data from RS Vul show a more feature-rich display, with out-of-plane velocities as high as $v_z = \pm 480$ km s$^{-1}$. These are much larger than the $\pm 210$ km s$^{-1}$ maximum we measure in our emissivity plots, but the magnitude of the tomogram velocities may be somewhat exaggerated due to the data's relatively low resolution in the z-direction. Another important difference from our model is that the intensities of the emission features in the tomograms remain on the order of those in the central plane even at very high values of $v_z$. The primary explanation for our lower intensity measurements is that the emission of low-density regions is suppressed by our use of an isothermal gas. Since the emissivity calculation depends strictly on density in the isothermal case, the magnitude of the emission in the simulation data appears smaller than it otherwise would in regions containing hot low-density gas. Since we have neglected the influence on the ionization state and temperature of the gas in our emissivity plots, we expect our intensities to fall drastically in the lower-density portion of the flow. In the early stages of the run, this results in the maximum intensities occurring at the location of the stream. The intensity elsewhere is smaller by a factor of $10^{-2}-10^{-4}$.

The late stages of our simulation display the accretion disc in a steady-state with a  morphology more akin to the tomographic data of U CrB in its disc mode \citep{agafonov09}. This tomographic data shows more concentrated emission in the Keplerian locus than the tomograms of RS Vul, as well as both a star-stream and stream-disc interaction region. Our late-stage emissivity plots are also dominated by emission from the Keplerian locus due to the presence of a disc, but the material driven away from the orbital plane in the stream-disc interaction remains visible. Due to the strong density dependence of the emissivity, the disc dominates the emissivity plots. Intensities elsewhere range from $10^{-3} - 10^{-7}$ below the maximum measurement on the central slice.

Our measurements of the tilt angle of the disc show that deflecting and boosting the stream introduces warping of up to $3^{\circ}$, but no significant tilt with respect to the orbital plane. Alternately, the lower-density gas blown upward above the orbital plane does display significantly larger tilt that measures up to $50^{\circ}$ at the outer edge of the disc. Obtaining a tilt measurement from the tomographic data for comparison is difficult, particularly for RS Vul, which does not display a well-defined accretion disc. In U CrB, z-velocities of $\pm 420$ km s$^{-1}$ were observed within the predicted locus of the Keplerian accretion disc. A tilted disc was inferred from the observation that emission in tomogram slices corresponding to negative z-velocities was stronger on the left side of the $v_y-v_x$ velocity-space plots, while positive z-velocity slices showed stronger emission in the lower-right side of the $v_y-v_x$ plot. Quantifying this tilt by assuming a z-velocity of $420$ km s$^{-1}$ and Keplerian velocities at the surface of the primary gives a tilt angle of $37.2^{\circ}$. Performing the same calculation at the outer Keplerian orbit produces a tilt of $56.1^{\circ}$, which is comparable to the tilt measured in our simulations for the plume of gas on the far side of the accretor. It should also be noted that the $v_z$ values are uncorrected for the lower resolution in the $v_z$ direction. \citet{agafonov06} claim the correction factor will be small, but any correction factor will serve to reduce the predicted tilt. Assuming a factor as large as $0.5$ would reduce the maximum z-velocities to $v_z = \pm 240$ km s$^{-1}$, which is close to the maximum $v_z$ we measure. This correction results in tilts of $20.8^{\circ}$ and $36.6^{\circ}$ at the inner and outer radii. Our simulations indicate that it is unlikely that the disc will be uniformly tilted, since the stream-disc interaction acts primarily on the outer edge of the disc. The tomographic data that depicts asymmetric intensities measured on opposite sides of the accreting star could be interpreted as warping in the disc, as opposed to a uniformly inclined disc. A large part of the out-of-plane flow in our simulations also occurs within the Keplerian locus, and could contribute to this measurement.

The primary difference between the deflected and undeflected stream in our simulations is that the modified stream enables large amounts of flow above and below the orbital plane. Boosting the stream allows gas to circle the primary without striking its surface, which in turn means the system retains angular momentum that would normally be destroyed in the stream-disc and star-stream collisions. The additional (non-z) angular momentum manifests in the out-of-plane flow, and is eventually dissipated in the accretion disc. The disc's dimensions appear larger and the densities slightly higher due to the absence of the star-stream interaction as an angular momentum dissipation mechanism. The bulk of the disc remains primarily within the orbital plane, although it is more dense and more uniform than the disc shown in the tomographic measurements of U CrB. We can reconcile the differences between the simulated and observed discs by noting that our isothermal approximation results in a gas that is more easily compressible. This will lead to a more compact accretion disc and the tendency for less gas to be blown away in the interaction regions. It is also possibile that the tomographic images of U CrB in a disc mode show the early stages of the disc, in which case the relative intensity of the disc would be lower, and the disc would have a more nonuniform structure.

We see no evidence of a transition between a stream and disc accretion mode in our simulations. This is likely a result of our use of a constant mass transfer rate through L1, and potentially a byproduct of our isothermal equation of state. Examining our emissivity plots at $2P_{orb}$ shows a morphology very similar to that at $18 P_{orb}$, but with the central slice intensity $1-2$ orders of magnitude higher than that on the remainder of the grid. While the disc is present, it is not fully formed, and appears less uniform than late in the run. A reduced rate of mass transfer could offer the disc material a chance to accrete without being replenished by material from the stream. A non-isothermal gas would contribute to this effect by resulting in a less compressible flow that would impede the formation of a stable disc.  During this time, we would expect to observe the system in its stream phase. Since the time necessary for a full disc to form  is large (preliminary runs showed $\dot{M}_{L1} = \dot{M}_{acc}$ was achieved after $40-50 P_{orb}$), variations in the mass transfer rate occurring across smaller timescales would prevent the full formation of a stable disc. Our use of an isothermal gas also increases the compressibility, and subsequently the density of our disc. Relaxing the isothermal approximation could result in a lower-density, less uniform disc.

\section{Conclusions}

Our simulations of RS Vul have shown that boosting the velocity of the accretion stream and deflecting it below the orbital plane at L1 introduces significant changes in the structure of the accretion flow. Using an isothermal equation of state, we simulate the formation of an accretion disc with both an undeflected stream and a stream boosted to Mach 30 and deflected by $45^{\circ}$. In the undeflected case, we measure the formation of an accretion disc that is symmetric about the orbital plane and displays less than one angular grid zone ($0.5^{\circ}$) of warping. Modifying the stream causes some radial enlargement of the disc, along with warping as large as $3^{\circ}$. In addition, a low-density plume of gas is ejected above the orbital plane on the far side of the accreting star with densities on the order of $10^{-2}\rho_{L1}-10^{-3}\rho_{L1}$ and tilt angles as high as $50^{\circ}$.

The intensities we measure in our pseudo-emissivity plots are strongest in the orbital plane due to our use of an isothermal gas, which causes the emissivity to depend solely on $\rho^2$. The increased out of plane flow that occurs due to the deflected stream increases the emissivity above and below the central plane by a minimum of five orders of magnitude. We anticipate that relaxing the isothermal condition would cause higher emissivities to be measured outside of the central plane, as in the original tomographic data.

Our results are consistent with the hypothesis that stream deflection is responsible for the ``tilted'' out-of-plane material that is measured in the tomograms of U CrB, and also for the extreme out-of-plane flow observed in RS Vul. We have shown that a strong velocity boost and deflection at L1 are sufficient to change the direct-impact nature of the system, which increases the amount of out-of-plane flow drastically. We expect that a variable mass transfer rate and potentially more extreme boosting would only serve to enhance the features observed in the out-of-plane flow.

\section{Acknowledgments}

This research was supported by an allocation of advanced computing resources provided by the National Science Foundation. The computations were performed on Kraken at the National Institute for Computational Sciences. This research was also supported in part by a North Carolina Space Grant award. The author would like to thank Dr. John Blondin for his assistance with this work, and the anonymous referee for the helpful suggestions regarding the manuscript.