\section{Introduction}

Solar energetic particles (SEPs) are electrons, protons and heavy ions that have been accelerated to supra-thermal energies (up to a hundred of MeV or even $\sim$~GeVs, in the case of ions), typically during solar eruptive events.
These charged particles may for example be energised at solar flare sites during magnetic reconnection events, or at shock waves driven by coronal mass ejections (CMEs).
After escaping from their acceleration site, SEPs propagate through the inner heliosphere, guided by the interplanetary magnetic field (IMF).
On their journey through space, a non-nominal solar wind as well as turbulent fluctuations may strongly affect the SEP trajectories, altering both the spatial and temporal dependencies of energetic particle distributions in the heliosphere.
Non-nominal solar wind conditions may for example be caused by interplanetary CMEs, or by solar wind streams of different speeds interacting with each other.
When a fast solar wind stream is followed by a slow stream, a rarefaction region forms between the two, whereas in the opposite case, the high-speed stream compresses the slow solar wind, which may eventually result in the formation of a corotating interaction region (CIR) \cite{smith76}.
CIRs may accelerate particles at large heliocentric radial distances, where one or both of the shock waves bounding the CIR are present, via diffusive shock acceleration \cite{bell78}, and at smaller radial distances due to the compressional waves bounding the transition region between the slow and the fast solar winds, e.g., \cite{desai98,giacalone02}.

To study the influence of a non-nominal IMF on SEP events, Wijsen et al.~\cite{wijsen19a} and \cite{wijsen19b} used a three-dimensional particle transport model to propagate particles in a solar wind generated by the EUropean Heliospheric FORecasting Information Asset (EUHFORIA) \cite{pomoell18}, which is a three-dimensional magnetohydrodynamic (MHD) model of the solar wind and CME evolution in the inner heliosphere.
Using EUHFORIA, these authors modelled a fast solar wind stream embedded in a slow solar wind, generating self-consistently a CIR from ${\sim} 1.3 $ AU onward.
By considering impulsive injections of 4 MeV protons from different extended injection regions located at 0.1 AU, Wijsen et al.~\cite{wijsen19b} illustrated how the longitudinal width of the particle intensity distributions can strongly increase or decrease as a function of radial distance, depending on the underlying IMF topology.
Moreover, by performing simulations both with and without cross-field diffusion, it was shown that  cross-field diffusion widens the region where particles propagate, i.e., the particle streaming zone, yet it does not influence the qualitative
behaviour of the longitudinal width as a function of radial distance \cite{wijsen19b}.
We want to point out that we leave it for future work to study the three-dimensional spread of particles in the early phase of SEP events, that is, where particles propagate closer to the Sun than our inner boundary at 0.1 AU.
Both observational and modelling studies indicate that particles may be accelerated by shocks associated with CMEs already $< 10$\,R$_{\odot}$ (e.g., \cite{roussev04,lee12,rouillard12,pomoell15,lario16,schwadron17,afanasiev18,kozarev19}, and references therein).
In addition, shock waves associated with CMEs contribute to the spread of SEPs in interplanetary space (e.g., \cite{cane06}).
Recently, Hu et al.~\cite{hu18} used the two-dimensional iPATH model (\cite{hu17}) to model both diffusive shock acceleration and interplanetary transport of protons starting from $10$\,R$_{\odot}$ during a gradual SEP event, as seen by virtual observers located at different radial distances in the solar equatorial plane, and covering a $60^\circ$ azimuthal range.
Assuming nominal solar wind conditions upstream of the CME driven shock wave, these authors used the Nonlinear Guiding Center (NLGC) theory \cite{matthaeus03,zank04,shalchi09} to  illustrate  how cross-field diffusion can alter the intensity-time profiles for 4.4~MeV and 46.4~MeV particles, especially for observers magnetically not well-connected with the travelling shock.
The intensity profiles measured by these observers showed more prompt enhancements for simulations with larger cross-field diffusion coefficients.

In this work, we extend the studies performed in \cite{wijsen19a} and \cite{wijsen19b} by applying different cross-field diffusion conditions to the particles propagating in non-nominal solar wind conditions in the interplanetary space.
In particular, we investigate how a larger mean free path perpendicular to the magnetic field might alter the results obtained in these previous studies.
After briefly introducing our energetic particle and solar wind models in  Section~\ref{sec:model}, we show how the longitudinal width of the particle intensity distribution varies as a function of radial distance, time, and level of cross-field diffusion (Section~\ref{sec:particle_spread}).
Next, we investigate how the different levels of cross-field diffusion used in our simulations affect the particle acceleration at the shock waves bounding the CIR (Section~\ref{sec:CD_acc}).
We end in Section~\ref{sec:discussion} with a summary and discussion of the obtained results.

\section{Simulation of SEP transport within a CIR}\label{sec:model}
To study SEPs in the inner heliosphere, we model the evolution of the directional particle intensity
$j(\vec{x},p,\mu,t)$
using the focused transport equation (FTE) (see e.g. \cite{roelof69, isenberg97,leRoux09}):
\begin{equation}\label{eq:fte}
\begin{aligned}
\pd{j}{t} &+\pd{}{\vec{x}}\cdot\left[\left(\td{\vec{x}}{t}+\pd{}{\vec{x}}\cdot \vec{\kappa}_\perp\right)j\right]+\pd{}{\mu}\left[\left(\td{\mu}{t}+ \pd{D_{\mu\mu}}{\mu}\right) j\right] + \pd{}{p}\left(\td{p}{t}j\right)\\ &= \pdd{}{\mu}\left[D_{\mu\mu}j\right] + \pd{}{\vec{x}}\cdot\left[\pd{}{\vec{x}}\cdot\left(\bm{\kappa}_\perp j\right)\right],
\end{aligned}
\end{equation}
with
\begin{eqnarray}
\td{\vec{x}}{t} &=& \vec{V}_{\rm sw}+ \mu \varv \vec{b} +\vec{V}_d  \label{eq:fte_x}  \\
\td{\mu}{t}&=&\frac{1-{\mu}^2}{2}\Bigg(\varv\nabla\cdot\vec{b}+ \mu \nabla\cdot\vec{V}_{\rm sw}- 3 {\mu} \vec{b}\vec{b}:\nabla\vec{V}_{\rm sw} - \frac{2}{{\varv}}\vec{b}\cdot\td{\vec{V}_{\rm sw}}{t}\Bigg)   \label{eq:fte_mu} \\
\td{p}{t} &=& \Bigg(\frac{1-3{\mu}^2}{2}(\vec{b}\vec{b}:\nabla\vec{V}_{\rm sw}) - \frac{1-{\mu}^2}{2}\nabla\cdot\vec{V}_{\rm sw} -\frac{{\mu} }{{\varv}}\vec{b}\cdot\td{\vec{V}_{\rm sw}}{t}\Bigg){p} \label{eq:fte_p}.
\end{eqnarray}
In these equations,
$\vec{x}$ denotes the phase-space spatial coordinate and $t$ the time, both measured in an inertial frame, whereas the cosine of the pitch angle $\mu$ and the momentum magnitude $p$ or speed $\varv$ are expressed in a frame co-moving with the solar wind.
Furthermore, $\vec{V}_{\rm sw}$ is the solar wind velocity, and
$\vec{b}$ the unit vector in the direction of the mean magnetic field,
$D_{\mu\mu}$ is the pitch-angle diffusion coefficient,
$\bm{\kappa}_\perp$ the spatial cross-field diffusion tensor, and
$\vec{V}_d$ the drift velocity due to the gradient and curvature of the mean magnetic field.
Our particle transport model solves Eq.~\eqref{eq:fte} by integrating the equivalent set of It\^o stochastic differential equations forward in time \cite{gardiner04}.
The solar wind velocity and magnetic field variables appearing in Eqs.~\eqref{eq:fte_x}--\eqref{eq:fte_p} are obtained from the three-dimensional MHD model EUHFORIA. For details regarding the numerical procedures used by our model, we refer to \cite{wijsen19a}.

As in \cite{wijsen19a} and \cite{wijsen19b}, we use EUHFORIA to model a slow solar wind with an embedded fast solar wind stream.
This is done by prescribing everywhere at the inner boundary of EUHFORIA a solar wind speed of 330~km/s, except for a region satisfying
 $
(\rm{longitude}-75^\circ)^2 + (\rm{latitude}-5^\circ)^2< (20^\circ)^2,
$
where we prescribe a solar wind of speed 660~km/s.
Moreover, this high-speed region is surrounded by a transition zone of width $6^\circ$.
This transition eventually evolves into a CIR, bounded by reverse and forward MHD shock waves.
As discussed in \cite{wijsen19a}, the precise location where the large-amplitude compression waves bounding the transition region evolve into shock waves, and hence marking the formation of the CIR, is difficult to pinpoint.
This is because the finite resolution of our MHD simulation produces inevitably shock waves that are wider than observed interplanetary CIR shocks.
However, we observe in our particle simulations that particles start to get accelerated from a heliographic radial distance of ${\sim}1.3$ AU onward, indicating that the large-amplitude compression waves have steepened sufficiently such that their characteristic length scale is considerably smaller than the particle mean free path that we assume (see below).

In contrast to the EUHFORIA simulation used in \cite{wijsen19a} and \cite{wijsen19b}, we assume in the current study a monopolar positive magnetic field at the inner boundary of the model, such that there is no current sheet present in our solar wind.
This is done to avoid any effects of the current sheet on the simulated particle intensities, which is necessary due to the larger spread of the particles in this work compared to the simulations performed in \cite{wijsen19a} and \cite{wijsen19b}.
We want to note that
the guiding centre approximation, used in Eq~\eqref{eq:fte}, breaks down near the current sheet \cite{speiser65}.
Apart from the magnetic field polarity, all other MHD variables prescribed at the inner boundary of the model are the same as in
\cite{wijsen19a}, and hence we refer to that paper for further details.

The simulated source of particles is located at the inner boundary of the model, at a heliographic radial distance of 0.1 AU, and centred around the equatorial plane with a latitudinal width of $10^\circ$.
The source has a longitudinal width of $30^\circ$ such that it covers a region extending from the slow to the fast solar wind.
As a consequence, the particle injection region is magnetically connected to both the forward and reverse compression/shock waves, hence also  encompassing the SI inside the CIR, where the shocked slow and fast solar wind plasma meet.
We note that the injection region considered in this work is identical to the one of \cite{wijsen19a} and to case~4 of \cite{wijsen19b} (see Fig~1 of \cite{wijsen19b}).
Likewise, we also consider a mono-energetic injection of 4~MeV protons, allowing us to easily track the energy changes of the particles during their propagation in the solar wind.
Unlike in \cite{wijsen19a} and \cite{wijsen19b}, in order to mimic a prompt injection of particles, we inject the protons according to a Reid-Axford time profile \cite{reid64}:
\[
j(t) \propto \frac{1}{t}\exp{\left(-\frac{\beta}{t} - \frac{t}{\tau}\right)},
\]
where we choose the parameters $\beta = 1.5 $~hours and $\tau=0.2$~hours (like in \cite{lario07}).

The particle scattering conditions are determined by the pitch-angle diffusion coefficient $D_{\mu\mu}$ and the cross field diffusion tensor $\bm{\kappa}_\perp$.
As detailed in  \cite{wijsen19a}, we use quasi-linear theory to prescribe the functional form of $D_{\mu\mu}$.
The strength of the magnetic turbulence acting as pitch-angle scattering centres is quantified by assuming a constant radial mean free path $\lambda_r^\parallel = 0.3$ AU for 4~MeV protons.
 For the cross-field diffusion tensor we assume an expression similar to the one used by e.g. \cite{zhang09,wang12}:
\begin{eqnarray}\label{eq:cross_field}
\bm{\kappa}_\perp &=&\alpha \frac{\pi }{12}{ {\varv}\lambda_\parallel}\frac{B_0}{B}\left(\mathbb{I}-\vec{bb}\right),
\end{eqnarray}
where  $\mathbb{I}$ is the unit tensor,
$\vec{bb}$ the dyadic product of the magnetic unit vector, and
$\lambda_\parallel$ the parallel mean free path.
$B_0$ is a reference magnetic field strength value fixed to  the maximum value of $B$ at 1 AU,
$B_0 = 9.7$~nT, which is attained in the shocked slow solar wind, adjacent to the SI.

We note that we used rather simple functional forms for the pitch-angle and cross-field diffusion coefficients.
However, the significant uncertainties regarding the turbulent conditions and hence the particle diffusion properties in large scale structures like CIRs, make it difficult to establish more realistic assumptions for the particle diffusion coefficients \cite{horbury99, crooker99,richardson18}.
Nonetheless, in Wijsen et al.~\cite{wijsen19a,wijsen19b} it is demonstrated how our choice for the pitch-angle and the cross-field diffusion coefficients can reproduce different features of the energetic particle intensity profiles associated with CIRs, like e.g. the particle acceleration near the shock waves and the intensity dip near the SI.
This is also discussed in more detail in Section~\ref{sec:discussion}.
We remark that in reality the parallel and perpendicular mean free path of the particles might vary rather strongly across the CIR due to the presence of different turbulent regimes \cite{horbury99}.
The increase in turbulence near the CIR shocks (see  e.g, \cite{tsurutani95} and Section~\ref{sec:discussion} for a discussion) would most likely decrease the parallel mean free path and amplify some of the effects that we describe in this work.
A smaller radial mean free path would e.g., imply that the particles will scatter across the CIR shock waves more often, resulting in more particle acceleration.
On the other hand, using a large constant radial mean free path everywhere in our simulations allows us to more easily trace the effects of the underlying non-nominal IMF topology and the different levels of cross-field diffusion.
In addition, we note that more advanced theories of the cross-field diffusion, like e.g., NLGC theory, require the specification of the magnetic turbulence level $\delta B$ in the heliosphere.
For nominal solar wind conditions, one typically assumes that this quantity decreases as a power-law with radial distance from the Sun (see e.g., \cite{mason12,zhao16, hu17}).
However, such a singular dependence on the radial coordinate is not valid in a solar wind containing a CIR, as the turbulent conditions in a CIR are significantly different from those encountered in the fast or slow solar wind (see e.g., \cite{horbury99, crooker99,richardson18}).
Therefore, since the behaviour of $\delta B$ in a CIR is not well understood, we instead adopt the phenomenological form~\eqref{eq:cross_field} for the cross-field diffusion, as its effects are more easily traceable in the simulations.

By defining the perpendicular mean free path as $\lambda_\perp = 3\kappa_\perp/\varv$, the parameter $\alpha$ determines the ratio $\lambda_\perp/\lambda_\parallel$ at 1 AU up to a factor of order unity.
We consider three different values for the parameter $\alpha$, namely $2\times10^{-4}$, $2\times10^{-3}$, and $2\times10^{-2}$. In the remainder of the text, we will refer to the simulations performed with these different values of $\alpha$ as case~1, case~2 and case~3, respectively.
These values of $\alpha$ help us to classify the `level' or `strength' of the cross-field diffusion employed in our simulations.
In this way, the case with the strongest cross-field diffusion corresponds to case~3 ($\alpha = 2\times10^{-2}$).

\section{Evolution of the particle spread in the solar equatorial plane}\label{sec:particle_spread}
\begin{figure*}
        \centering
        \begin{tabular}{cc}
        \includegraphics[width=0.45\textwidth]{dlon_evolution_1K}
        &
        \includegraphics[width=0.45\textwidth]{lonBoundaries_1K}\\
        \includegraphics[width=0.45\textwidth]{dlon_evolution_10K}&
        \includegraphics[width=0.45\textwidth]{lonBoundaries_10K}\\
        \includegraphics[width=0.45\textwidth]{dlon_evolution_100K}&
        \includegraphics[width=0.45\textwidth]{lonBoundaries_100K}\\
    \end{tabular}
    \caption{Column~A (left column) shows the longitudinal width of the omni-directional intensity distribution as a function of radial distance for the different indicated times.
    Column~B (right column) shows a snapshot at 11.5~hours of the omni-directional particle intensities, drawn in grey shades on top of the radial velocity profile of the solar wind.
    The green and blue curves represent, respectively, the lower and upper longitudinal boundaries of the streaming zone, as determined by the prescribed intensity threshold value $10^{-8.5}$ (arbitrary units).
 The first row is for case~1, the second row for case~2, and the third row for case~3. }
     \label{fig:dlon}
\end{figure*}
In this section, we analyse how the different levels of cross-field diffusion affect the longitudinal width of the particle intensity distributions in the solar equatorial plane.
For this purpose, we define the particle streaming zone as the region where the omni-directional intensity, integrated over energy, is larger than the limiting value $10^{-8.5}$ (in arbitrary units).
This threshold is set three orders of magnitude below the peak-intensity attained at 1 AU for a simulation without cross-field diffusion.
Defining the longitudinal coordinate $\varphi$ in anticlockwise
direction, we obtain at each heliographic radial distance, $r$, and time, $t$, the lower, $\varphi_L(t,r)$, and upper, $\varphi_U(t,r)$, edge longitudes of the particle streaming zone, such that we can define the longitudinal width of the particle streaming zone as $\Delta\varphi(r,t) =\varphi_U(t,r) - \varphi_L(t,r)$.
The left column of  Fig.~\ref{fig:dlon} shows for the three cases the radial dependence of $\Delta\varphi(r,t)$ at different times.
The right column of Fig.~\ref{fig:dlon} shows for each case the corresponding snapshots of the particle intensities at 11.5~hours.
In these figures, the lower longitudinal boundary is indicated by a green line, whereas the upper longitudinal boundary is  shown in blue.

From the first column of Fig.~\ref{fig:dlon}, we see that for all cases, $\Delta\varphi(r,t)$ generally decreases with increasing radial distance, reflecting the IMF compression at the shock waves bounding the CIR (see also Fig.~6 of \cite{wijsen19a}).
However, apart from this common decreasing trend, there are several differences between the radial evolution of $\Delta\phi$ when comparing the different cases.

For case~1, panel~1A shows that at large radial distances ($r\gtrsim 2$ AU), the width of the particle streaming zone remains approximately constant and even starts to increase slightly.
At these distances, the particle streaming zone is located entirely within the CIR, meaning that all the IMF lines on which particles are travelling have crossed the shock waves at distances $r<2$ AU.
In contrast, we see in panel~2A that for case~2 the longitudinal width of the particle streaming zone does not become constant from $r=2$ AU onward.
This is because the stronger cross-field diffusion has brought particles to IMF lines that cross the shock waves at radial distances larger than $2$ AU, and hence, contrary to case~1, there are still particles outside the CIR at these large radial distances.
However, since the forward and reverse shock waves propagate into the slow and fast solar wind, respectively, we expect that from a certain radial distance onward, all particles will again reside inside the CIR, and hence the behaviour of the longitudinal width will probably become qualitatively similar to what is observed for case~1.

Panel~3A shows that the strong cross-field diffusion of case~3 introduces several features to the radial dependence of the longitudinal width that are not seen in the two other cases.
At small radial distances, the width of the particle streaming zone remains approximately constant, and even increases slightly, which is because some particles have been transported to the unperturbed slow solar wind behind the high-speed stream.
However, panel 3A shows that the constant longitudinal width of the particle intensity distribution is followed by a sudden drop at larger radial distances.
For times later than 9.5 hours, this drop occurs at approximately the same radial distance, $r=1.4$~AU.
Looking at panel 3B, we see that this decrease is because the omni-directional intensity drops below the $10^{-8.5}$ threshold in the rarefaction region behind the high-speed solar wind stream (the pronounced indent in the green intensity threshold contour in Figure 3B).
This is a result of the strong expansion of magnetic flux tubes inside the rarefaction region, which causes the energetic particle density to get rarefied as well.
In addition, the large cross-field diffusion for case~3 spreads the energetic particles over a much larger region than for the two other cases, hence rarefying the particle densities everywhere in the heliosphere.
This is especially true in the fast solar wind stream and the rarefaction region, since the low magnetic field strength implies a more efficient cross-field diffusion (see Eq.~\eqref{eq:cross_field}).
As a consequence, the intensities for case~3 are everywhere relatively low as compared to the previous cases.
These low intensities might imply that a large part of the SEP event may go undetected if  the pre-event intensity levels are high enough or the measuring instrument is not sensitive enough.
In fact, such a situation is illustrated in panel~3B, where we see that a large region is characterised by particle intensities below the $10^{-8.5}$ threshold, which is therefore not taken into account when determining the width of the intensity distributions.
The low intensity values due to the efficient particle spreading also explains why in panel~3A the longitudinal width does not increase any further after 9.5~hours, and why for a radial distance larger than 1.5~AU, the longitudinal width at time 15.5~hours is even smaller than at the three previous time instances.

From the first column of Fig.~\ref{fig:dlon}, we note that for all three cases the longitudinal width of the particle streaming zone at small radial distances and at 1.5~hours is already significantly larger than the longitudinal width of the injection region, which was $30^\circ$.
A rapid spread of particles at early times can be understood by considering the diffusive cross-field flux, which can, for  small radial distances and in the solar equatorial plane, be approximated as in \cite{strauss17b}:
\[
\left|\mathcal{F}_\perp^{\rm dif}\right|\sim \kappa_\perp  \left|\pd{j}{\varphi} \right|.
\]
From this expression, we see that the cross-field diffusion will be more efficient when large longitudinal gradients are present in the directional intensity function.
Since our particles are injected with a uniform intensity from a source region of limited spatial extend, we have in fact that, upon injection,  $\left|\partial{j}/\partial{\varphi} \right| \longrightarrow +\infty$ at the boundaries of our injection region in the solar equatorial plane.
We note that this infinite derivative does not produce numerical problems due to the stochastic nature of our model.
However, it explains the efficient cross-field diffusion early on in the event (see also \cite{strauss17b} for a discussion).

Moreover, we note that the above argument is more generally valid at any location in the heliosphere when considering the gradient of the intensity in the direction perpendicular to the magnetic field.
More specifically, for an isotropic $\bm{\kappa}_\perp$, the diffusive cross-field flux in the direction of a unit vector $\vec{e}_\perp$, with $\vec{e}_\perp \perp \vec{b}$, can be defined as
\begin{equation}\label{eq:cross-fieldFlux}
  \left|\mathcal{F}_\perp^{\rm dif}\right|\sim \kappa_\perp  \left|\vec{e}_\perp\cdot\nabla{j}\right|.
\end{equation}
Hence, cross-field diffusion will be most effective where $\kappa_{\perp}$ is large or where there are large particle intensity gradients perpendicular to the magnetic field.
Such large gradients are induced by the IMF topology inside the CIR, where the magnetic field is compressed.
As discussed above, this compression of the magnetic field lines leads to a compression of the particle intensity distribution, and hence a decreasing longitudinal width.
However, the same compression will also produce larger particle intensity gradients perpendicular to the magnetic field, hence increasing the effectiveness of cross-field diffusion according to Eq. \eqref{eq:cross-fieldFlux}.
At the same time, we note that the inverse scaling of our cross-field diffusion with the magnetic field strength implies that the perpendicular mean free path decreases inside the CIR.
However, as discussed in \cite{wijsen19a} and \cite{wijsen19b}, a small cross-field motion of the particles inside the CIR where the IMF lines are closely compressed may transport the particles on field lines that are significantly separated outside the  CIR, i.e., at smaller radial distances.
The observed cross-field diffusion in our simulations results thus from a non-trivial interplay between the magnetic field dependence of the cross-field diffusion tensor and the intensity gradients perpendicular to the magnetic field.
Both mechanisms are strongly influenced by non-nominal solar wind conditions.

Finally, we remark that by adding the cross-field diffusion term in the spatial part of the FTE, i.e., in Eq.~\eqref{eq:fte_x}, the energetic particles  are given an additional velocity component on top of their actual velocity $\vec{\varv}$.
 As a consequence, particles with $\mu \sim 1$ might spread faster than their speed allows for, hence violating causality  \cite{strauss15}.
 Following a similar procedure as in~\cite{strauss15} to estimate the impact of this effect in our simulations, we found it to be negligible since the intensities resulting from such particles are more than two orders of magnitude below the selected threshold value.
This non-physical behaviour is a consequence of treating the motion of particles perpendicular to the mean magnetic field as a spatial diffusion process, and would not appear if one would instead e.g., model cross-field motions by stochastically varying only the propagation direction of the particles.
Such an approach boils down to  varying the magnetic unit vector, which is similar to the method implemented by Laitinen et al.~\cite{laitinen13} to model the  motion of particles along random-walking magnetic field lines.

\section{The effect of cross-field diffusion on shock accelerated particle populations}\label{sec:CD_acc}

\begin{figure*}
        \centering
        \begin{tabular}{cc}
        \includegraphics[width=0.45\textwidth]{protons_CIR_noK_R1_00auCLT90_00deg_T11_000dn0_60_S0_75_1K}&
        \includegraphics[width=0.45\textwidth]{protons_CIR_noK_R1_00auCLT90_00deg_T31_000dn0_60_S0_75_1K}
    \end{tabular}
    \caption{Snapshot of the proton intensities (upper panels) and anisotropies (lower panels) for case~1  at 11.0~hours (left panels) and 31.0~hours (right panels). The profiles are measured in the solar equatorial plane at a fixed heliocentric radial distance of 1 AU, and shown as a function of the longitudinal separation from the centre of the injection region.}
     \label{fig:lon1}
\end{figure*}

\begin{figure*}
        \centering
        \begin{tabular}{cc}
        \includegraphics[width=0.45\textwidth]{protons_CIR_noK_R1_00auCLT90_00deg_T11_000dn0_60_S0_75_10K}&
        \includegraphics[width=0.45\textwidth]{protons_CIR_noK_R1_00auCLT90_00deg_T11_000dn0_60_S0_75_100K}\\
    \end{tabular}
    \caption{Snapshot of the proton intensities (upper panels) and anisotropies (lower panels) at time 11.0~hours  for case~2  (left panel) and case~3  (right panel). The profiles are measured in the solar equatorial plane at a fixed heliocentric radial distance of 1 AU, and shown as a function of the longitudinal separation from the centre of the injection region.}
     \label{fig:lon2}
\end{figure*}

We now examine the particle intensity and anisotropy profiles in the solar equatorial plane at a fixed heliographic radial distance of 1~AU.
We focus hereby on energy channels above the initial 4~MeV injection energy of the particles, meaning that we only consider particles that have been accelerated at the shock waves bounding the CIR.
We note that the acceleration of the protons and ions in shock waves depends on several factors, including the scattering conditions (e.g., \cite{lee83,gordon99, vainio07,afanasiev15}), the shock strength (e.g., \cite{lario98,rice03,gasen11,gasen14,prinsloo19}), and the shock obliquity (e.g., \cite{zank01,zank06,battarbee13}), among others.
In Section~\ref{sec:discussion} we provide a discussion on these factors in the context of CIR shocks, and how they influence the results presented in this section.

Figure~\ref{fig:lon1} shows two snapshots, at 11.0~hours (left) and at 31.0~hours (right), of the intensity and anisotropy profiles for case~1, as a function of $\varphi - \varphi_c$, where $\varphi_c$ denotes the longitude of the point at $1$ AU that connects magnetically to the centre of the particle injection region at 0.1 AU.
The left panel of Fig.~\ref{fig:lon1} shows two distinct intensity enhancements in the $4.1\pm 0.1$ MeV energy channel, resulting from particle acceleration at the two shock waves bounding the CIR.
The intensity enhancement at $\varphi - \varphi_c <0$ is produced by particle acceleration at the reverse shock, whereas the intensity enhancement at $\varphi - \varphi_c >0$ is produced by particle acceleration at the forward shock.
The negative anisotropies in the lower panels of Fig.~\ref{fig:lon1} indicate that the acceleration of the particles occurs at radial distances larger than 1~AU.
This is in accordance with the findings of Wijsen et al.~\cite{wijsen19a}, where we observed the generation of accelerated particles at the shocks for $r \gtrsim$ 1.3 AU.
In our simulations, the reverse shock is a considerably more efficient particle accelerator than the forward shock.
For the former, we see that after 11~hours there are particles accelerated up to ${\sim}1.0$~MeV above their initial energy, whereas the latter only accelerates particles slightly above 4~MeV.

The two intensity enhancements in the left panel of Fig.~\ref{fig:lon1} can co-exist without merging partly because the SI separating them is characterised by a strong magnetic field, and hence by a weak cross-field diffusion since the latter scales inversely proportional to the magnetic field strength (see Eq.\eqref{eq:cross_field}).
However, since the cross-field diffusion is not exactly zero at the SI, the separation between the two intensity enhancements will slowly diffuse away.
The latter might be aided by the large particle intensity gradients perpendicular to the magnetic field observed at the SI.
Looking at the right panel of Fig.~\ref{fig:lon1}, which shows the intensity and anisotropy profiles of case~1 at 31.0~hours, we realise that the intensity dip has indeed increased to a higher intensity level. By this time, the reverse shock has accelerated protons above 5.3~MeV, although the intensities and hence the statistics for these particles are low.

Figure.~\ref{fig:lon2} shows the longitudinal intensity and anisotropy profiles at 1 AU for case~2 and~3 at 11.0~hours.
Focusing on the left  panel of Fig.~\ref{fig:lon2}, we see that the stronger cross-field diffusion of case~2 has merged the two $4.1\pm 0.1$ MeV intensity enhancements of case~1 into one wide enhancement.
Moreover, we see that there is less particle acceleration when compared to the previous case, since the
$5.1\pm 0.1$ MeV and $5.3\pm 0.1$ MeV energy channels are no longer populated by particles.
The peak intensities of the other channels have also decreased, yet we note that these channels cover now a larger longitudinal range.
The drop in accelerated particle intensities with increasing cross-field diffusion is even more notable when looking at the right panel of Fig.~\ref{fig:lon2}, which shows the intensity profiles for case~3.
This panel illustrates that the shock waves are barely capable of accelerating particles due to the large cross-field diffusion assumed in case~3.
This is because the cross-field diffusion efficiently moves particles away from the shock waves by transporting the particles to IMF lines that either cross the shock wave at much larger radial distances, or that enter the CIR before the shock waves have formed.

\section{Summary and discussion}\label{sec:discussion}

In this work, we have extended the studies performed in \cite{wijsen19a} and \cite{wijsen19b} that investigate SEP transport near and within a high-speed solar wind stream that eventually evolves into a CIR.
For this, we used a three dimensional particle transport model that computes particle directional intensities and anisotropies in a solar wind generated by EUHFORIA.
In particular, we considered an injection of 4~MeV protons at the inner boundary of our model, from a source region that extends both in the slow and fast solar wind regimes, hence also covering the transition region between both plasma environments.
Assuming a cross-field diffusion coefficient that scales inversely proportional to the magnetic field strength, we performed simulations with three different levels of cross-field diffusion.
This allowed us to study how the strength of cross-field diffusion may alter the particles spread and the energy content of SEPs in a solar wind that deviates from a nominal solar wind configuration, under the assumptions of our model.

From our simulations, we found that the longitudinal widths of the particle intensity distributions show a general decreasing trend with increasing radial distance, reflecting the underlying IMF topology near and within the CIR.
Such a decreasing width is in contrast with what is observed in nominal solar wind configurations, where the width of the particle intensity distribution remains constant in the absence of cross-field diffusion or increases with radial distance if cross-field diffusion is present (see \cite{wijsen19b}).
Although our three simulations with different levels of cross-field diffusion showed overall a similar decreasing trend, there were some notable differences.
For the case with the strongest cross-field diffusion, there was a sudden drop in the longitudinal width of the particle intensity distributions at ${\sim}1.4$~AU, because the particle intensity falls below our prescribed threshold.
This intensity drop occurred in the rarefaction region behind the high-speed solar wind stream, where the magnetic flux tubes are strongly expanding, hence rarefying the SEP densities.
In addition, the strong cross-field diffusion in our simulation dispersed particles over such a large region in the heliosphere that the intensities decreased everywhere significantly in comparison with the simulations that use a weaker cross-field diffusion.
As a consequence, a substantial amount of particles resided in regions of the solar wind showing intensities below our prescribed threshold, and were therefore not taken into account when calculating the longitudinal width.
A similar situation may happen with actual in-situ observed SEP events, i.e., when particles propagate either under strong diffusive transport conditions or in solar wind rarefaction regions, they might be dispersed efficiently such that their intensities drop below the background level measured by a given particle instrument, as similarly noticed by Pacheco et al.~\cite{pacheco19} for electron events simulated under diffusive parallel transport conditions.

Similar to what Strauss et al.~\cite{strauss17b} obtained for electron propagation in a nominal IMF, we find that the cross-field diffusion is very effective early in the event, when large intensity gradients are present.
However, in contrast to a Parker solar wind, the compressed IMF topology inside the CIR also gives rise to large intensity gradients perpendicular to the magnetic field, increasing the effectiveness of cross-field diffusion.
In contrast, the strongly diverging magnetic field topology in the rarefaction region might produce very low particle intensity gradients, hence decreasing the diffusive cross-field flux as defined in Eq.~\eqref{eq:cross-fieldFlux}.
We find thus that for a cross-field diffusion coefficient that scales inversely proportional to the magnetic field strength, the effectiveness of the cross-field diffusion is determined by an intricate interplay between the variations in the magnetic field strength and the particle intensity gradients, both modulated by the IMF topology and hence by the converging or diverging properties of the solar wind flow.

In Section~\ref{sec:CD_acc}, we showed the intensity profiles of the shock-accelerated particles at 1~AU, as a function of longitude.
For case~1, i.e., the case with the weakest cross-field diffusion, two separate intensity enhancements were visible, reflecting the particle acceleration occurring at the two different shock waves.
The negative anisotropies indicate that the acceleration of the particles happens at larger radial distances, where the shocks have formed.
This agrees with in-situ observations at 1 AU, showing that the flow of energetic ions accelerated at a CIR shock is typically towards the Sun \cite{marshall78}.

Moreover, the double structure in the intensity profiles is a feature that is also commonly detected with in-situ observations at different radial distances from the Sun \cite{barnes76,decker81,tsurutani82,dwyer97,mason99}.
Like in our simulations, observations measure typically more accelerated particle fluxes originating from the reverse shock than from the forward shock \cite{tsurutani82,desai98,mason99b}.
This difference might be partly due to the fact that the relative speed difference between the unperturbed solar wind and the shock is higher for the reverse shock than for the forward shock, since the reverse shock propagates radially inwards in the solar wind frame (see also \cite{03kocharov}).

Another important factor that might contribute to the observed differences in  accelerated particle fluxes originating from the different shocks is the angle between the shock normal and the upstream magnetic field, i.e. the shock angle, $\theta_{Bn}$.
In our simulation, the reverse shock of the CIR is more oblique than the forward shock, which follows from the difference in speed between the fast and slow solar wind.
Tsurutani et al.~\cite{tsurutani82} demonstrated that  proton fluxes near CIRs
are correlated with the shock angle, with the most intense proton events occurring for shocks that are quasi-perpendicular.
This is in accordance with the results obtained by Zank et al.~\cite{zank06}, which illustrated that perpendicular shock can be very efficient particle accelerators given that there is a sufficient amount of turbulence present in the solar wind to diffuse the particle across magnetic field lines.
Interesting to note is that CIR shocks are indeed typically characterised by increased amounts of turbulence (see e.g. \cite{horbury99}) as compared to the unperturbed solar wind, despite the fact that wave excitation by the streaming instability is quenched for highly oblique shocks.
The latter is a consequence of the proportionality of the  wave growth to the shock angle \cite{gordon99,li03,rice03}.

The dip between the two intensity enhancements at the reverse and forward shocks coincides with the SI,
and can be explained by noting that the field lines adjacent to the SI do not cross the reverse shock wave, nor the forward shock wave.
Hence the only way for shock-accelerated particles to reach these field-lines is through cross-field transport.
Using observations made by the Pioneer~10 and~11 spacecraft of a CIR lasting for 14 rotations, Tsurutani et al.~\cite{tsurutani82} showed that the intensity dip increases with a factor of 10 over 11 solar rotations, likely due to cross-field diffusion.
Similarly, we see in the simulations of case~1 that the intensity minimum increases in time.
Figure~\ref{fig:lon1} illustrated that the dip increased by a similar factor of ${\sim}10$, yet in only 20 hours, i.e., much faster than the observations in \cite{tsurutani82}.
This difference can have various causes, as discussed in the following.
It may suggest that the cross-field diffusion used in case~1, for which $\lambda_\perp/\lambda_\parallel{\sim} 10^{-4}$ at 1 AU, might still be too strong inside a CIR.
However, Dwyer et al.~\cite{dwyer97} obtained values for $\lambda_\perp/\lambda_\parallel$ of the order of unity for three CIRs using 44\,--\,313~keV/nucleon helium measurements from the Wind spacecraft.
Important to note here is that the results of Tsurutani et al. \cite{tsurutani82} regarding the increase in the intensity dip near the SI were done when the Pioneer~10 and~11 spacecraft were at a radial distance of ${\sim} 4.4$ and $4.9$ AU, whereas  the Wind spacecraft is located at ${\sim} 1 $ AU, and hence the properties of the solar wind turbulence might be different.
Other non-negligible factors that affect the acceleration of particles and hence also the evolution of the intensity dip are the plasma normal velocity jump and the compression ratio across the shocks as well as the obliquity of the shocks.
In fact, Desai et al.~\cite{desai98} showed that the proton intensity at CIRs is well correlated with the compression ratio, even if there is no shock present, at least when the particle seed remains constant.
In our simulations the particles accelerated at the CIR shocks originate from a flare-like impulsive event, with the particle source region magnetically connected to a limited part of the CIR shocks.
However, studies have shown that the seed particle population of most CIR events is likely the supra-thermal tail of the solar wind (see e.g. \cite{chotoo00,mason12b}).
This means that the seed population consists of charged particles with energies around $\sim 10$ keV/nucleon, i.e., much lower than the 4 MeV injection energy that we consider.
Pitch-angle scattering, cross-field diffusion and the acceleration processes in these non-nominal solar wind conditions are dependent of the momentum of the particles.
A continuous reservoir of low-energy seed particles at the CIR shocks may thus also help explaining the differences observed between our model and the observations in \cite{tsurutani82}.
However, the observation that CIRs accompanied by the largest energetic particle enhancements are often associated with quasi-perpendicular shocks~\cite{tsurutani82}, raises the question how low energy protons would overcome the associated injection problem \cite{jokipii87,zank96}.
Possible explanations can be found in the high levels of turbulence found near CIR shocks, or by considering a seed population originating from the supra-thermal tail  of a shock-heated solar wind plasma.
Prinsloo et al.~\cite{prinsloo19} illustrated that the  prior heating of the solar wind  distribution is important in reproducing the intensities of typical energetic storm particle events observed at shock waves driven by CMEs, especially when the injection energy is high like is typically the case for quasi-perpendicular shocks \cite{tylka05,zank06}. A similar conclusion may hold for the shock waves bounding the CIR.
We note that in our simulations we avoid any injection problem by considering protons of sufficiently high energy only.
We leave it for future work to study how our results change when e.g., continuously injecting the low-energy particles from the supra-thermal tail of the solar wind.

The double-enhancement intensity structure never forms at 1 AU in our case~2 and~3 simulations, where we use a cross-field diffusion magnitude comparable to that of Dwyer et al.~\cite{dwyer97}. This may be due to strong turbulence variations near and within the CIR that are not captured by our cross-field diffusion model.
More specifically, the plasma near the reverse shock is often characterised by significant magnetic fluctuations \cite{tsurutani95, horbury99}, while inside the CIR, the turbulent fluctuation of the magnetic field typically decreases towards the SI \cite{tsurutani82,dwyer97}.
More recently, Strauss et al.~\cite{strauss16} showed the decrease in turbulence at the SI by computing the variance of the in-situ measured magnetic field components for a CIR.
In addition, these authors argue that a tangential magnetic discontinuity like the SI may produce an anisotropic cross-field diffusion tensor, by strongly damping magnetic field fluctuations in the direction perpendicular to the discontinuity.
As a consequence, the SI would act as an efficient diffusion barrier, while still allowing for strong cross-field diffusion near the CIR boundaries.
Hence, the model proposed by Strauss et al.~\cite{strauss16} seems to explain both the slow increase of the intensity dip in the CIR shown in \cite{tsurutani82} and the large values of $\lambda_\perp/\lambda_\parallel$ reported in \cite{dwyer97} near the reverse shock.
In our simulations, a strong cross-field diffusion near the shocks leads to less particle acceleration, owing to the efficiency of the cross-field diffusion at spreading particles away from the shock wave.

As mentioned above, in addition to different shock properties, this apparent inconsistency between our simulation results and in-situ observations may be due to the fact that we are considering 4 MeV protons, whereas Dwyer et al.~\cite{dwyer97} obtain their results for helium with lower energies per nucleon.
The ion parallel mean free path typically decreases with decreasing rigidity \cite{droge00,droge03}, and therefore low-energy particles would be more efficiently trapped near the shock waves.
In addition, the strong levels of turbulence observed near CIR shocks might reduce the parallel mean free path of the particles significantly, hence both increasing $\lambda_\perp/\lambda_\parallel$ and facilitating  multiple shock encounters.
Moreover, we note that cross-field diffusion might have a non-negligible pitch-angle dependence, as discussed in \cite{strauss15}, which could also influence the scattering of the particles near the shock waves.
Further investigation of these issues is required, as it will eventually lead to a better understanding of particle transport in the complex turbulent environments that characterise large scale interplanetary structures like CIRs. The strongly changing turbulence properties found in the different parts of a CIR, and the possibility to have in-situ measurements of these regions, make CIRs the ideal laboratories to study the transport of particles perpendicular to the mean magnetic field.

\section*{Acknowledgments}
N. Wijsen is supported by a PhD Fellowship of the Research Foundation Flanders (FWO). The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders (FWO) and the Flemish Government, department EWI. The work at KU Leuven was done in the framework of the projects GOA/2015-014 (KU Leuven), G.0A23.16N (FWO-Vlaanderen) and C 90347 (Prodex). Funding for this work at University of Barcelona was partially provided by the Spanish MINECO under project MDM-2014-0369 of ICCUB (Unidad de Excelencia `Mar\'ia de Maeztu'). The work at University of Helsinki was carried out in the Finnish Centre of Excellence in Research of Sustainable Space (Academy of Finland grant numbers 312390 and 312351).

\section*{References}
\input{main.bbl}

\end{document}


\section{Introduction}

Stochastic reachability is an established verification tool to assure that a system will reach a desired state without violating predefined safety constraints
with at least a desired likelihood.
The solution to the stochastic reachability problem is based on dynamic programming
\cite{summers2010verification, abate2008probabilistic},
which poses significant computational challenges.
Optimization-based solutions have garnered modest computational tractability
via chance constraints \cite{lesser2013stochastic, vinod2019piecewise}, sampling
methods \cite{sartipizadeh2018voronoi, vinod2018multiple, vinod2018stochastic}, and convex optimization with Fourier
transforms \cite{vinod2018scalable, vinod2017scalable}, but are limited to linear
dynamical systems with Gaussian or log-concave disturbances.
Methods using approximate dynamic programming
\cite{kariotoglou2013approximate} and
particle filtering \cite{manganini2015policy, lesser2013stochastic},
are applicable to classes of nonlinear systems, but have only been demonstrated on systems with
moderate dimensionality.

Further, for many dynamical systems, presumption of accurate knowledge of dynamics and uncertainty is unrealistic. Historically, such uncertainty is handled through approximations and introduction of error terms that bound unknown elements \cite{mitchell2005time, fisac2015reach}.  With the rapid increase in the use of learning elements that are resistant to traditional models for control and formal methods, as well as the involvement of humans, such approaches may either be overly conservative or even simply inaccurate. For example, human inputs may be highly heterogeneous, may not follow a known distribution, and are often data-driven processes that may be biased when analyzed through sampling methods.

\emph{We propose to use conditional distribution embeddings within a reproducing kernel Hilbert space
(RKHS) to solve the stochastic reachability problem.}
As a tool for stochastic reachability, kernel methods offer significant advantages over state-of-the-art: 1) they are model-free, so can accommodate data-driven stochastic processes and nonlinear dynamics, 2) they are convergent in probability, and 3) they do not suffer from the curse of dimensionality
\cite{bellman2015applied}
or from issues of numerical
quadrature that plague optimization based approaches.
The primary computational challenge arises in the inversion of a matrix that scales with the number of data points, leading to computational complexity that is exponential in the size of the data.

Kernel methods are an established learning technique
\cite{scholkopf2001learning, shawe2004kernel, aronszajn1950theory}, and have
recently been applied to dynamic programming problems
with additive cost function and infinite time horizon
\cite{grunewalder2012modelling}.
Kernel methods broadly enable nonparametric inference using
kernel embeddings of distributions.  They can
capture the features of arbitrary statistical
distributions in a data-driven fashion \cite{song2013kernel,
smola2007hilbert}.
Kernel methods
do not suffer from biases or prior assumptions on the system model, and
are computationally efficient because they are convergent and
non-iterative \cite{scholkopf2001learning}.
These techniques have also been applied to several problems in dynamical
systems, including controller synthesis
\cite{lever2015modelling}, partially-observable systems \cite{nishiyama2012hilbert}, and estimation
of graphical models \cite{song2009hilbert}.

{\em The main contribution of this paper is the application of conditional distribution embedding
to compute the stochastic reachability probability measure, to enable model-free verification without
invoking a statistical approach.}
This is particularly relevant for systems with
black-box elements, such as autonomous or human-in-the-loop systems, which have
traditionally been resistant to formal verification techniques.
We tailor the approach in \cite{grunewalder2012modelling} to
accommodate the multiplicative cost and finite horizon associated with the stochastic reachability
dynamic program.
We apply kernel methods to compute the stochastic reachability
probability measure
by representing
the stochastic kernel
as a conditional distribution embedding within an RKHS.

The paper organization is as follows.
Section \ref{section: problem formulation} formulates the problem.
Section \ref{section: kernel mean embeddings} applies
conditional distribution embeddings to compute the stochastic reachability probability measure.
Section \ref{section: numerical results} demonstrates our approach on
three examples: a double integrator to enable validation with a ``truth'' model
via dynamic programming, a 10,000-dimensional integrator to demonstrate scalability, and spacecraft rendezvous and docking.

\section{Problem Formulation}
\label{section: problem formulation}

For sets $\mathcal{A}$ and $\mathcal{B}$,
the set of all elements of $\mathcal{A}$ which are not in $\mathcal{B}$
is denoted as $\mathcal{A} \backslash \mathcal{B}$.
We denote the indicator function
as $\boldsymbol{1}_{\mathcal{A}}(x) = 1$ if $x \in \mathcal{A}$
and $\boldsymbol{1}_{\mathcal{A}}(x) = 0$ if $x \notin \mathcal{A}$.

Let $\Omega$ denote a sample space
and $\mathcal{F}(\Omega)$ denote the $\sigma$-algebra relative to $\Omega$.
A probability measure $\Pr$ assigned to
the measurable space $(\Omega, \mathcal{F}(\Omega))$
is defined as the probability space $(\Omega, \mathcal{F}(\Omega), \Pr)$.
When $\Omega \equiv \Re$,
the $\sigma$-algebra of $\Omega$ is denoted as $\mathscr{B}(\Omega)$,
and is the Borel $\sigma$-algebra associated with $\Omega$.
A random variable $\boldsymbol{x}$ is a measurable function
on the probability space $(\Omega, \mathcal{F}(\Omega), \Pr)$.
A random vector $\boldsymbol{x} = [\boldsymbol{x}_{1}, \ldots, \boldsymbol{x}_{n}]^{\top}$
of $n$ random variables $\lbrace \boldsymbol{x}_{i} \rbrace_{i=1}^{n}$, each measurable functions on the probability space $(\Omega, \mathcal{F}(\Omega), \Pr)$, is defined on the induced probability space
$(\Omega^{n}, \mathcal{F}(\Omega^{n}), \Pr_{\boldsymbol{x}})$,
where $\Pr_{\boldsymbol{x}}$ is the induced probability measure.
A stochastic process is defined as a sequence of random vectors
$\lbrace \boldsymbol{x}_{k} : k \in [0, N]\rbrace$, $N \in \mathbb{N}$,
where $\boldsymbol{x}_{k}$ are defined on the probability space
$(\Omega^{n}, \mathcal{F}(\Omega^{n}), \Pr_{\boldsymbol{x}})$.
See \cite{billingsley2008probability, chow2012probability} for more details.

The expectation operator is denoted as $\mathbb{E}[\,\cdot\,]$,
where for some function $f$,
$\mathbb{E}_{\boldsymbol{x} \sim \Pr_{\boldsymbol{x}}\lbrace\,\cdot\,\rbrace}[f(\boldsymbol{x})]$
denotes the expectation operator
with respect to the probability measure $\Pr_{\boldsymbol{x}}$.

\subsection{Terminal-Hitting Time Problem}
\label{section: system model}

Consider a Markov control process $\mathcal{H}$,
which is defined in \cite{summers2010verification} as a 3-tuple,
\begin{align}
	\mathcal{H} = ( \mathcal{X}, \mathcal{U}, Q )
	\label{eqn: Markov control process}
\end{align}
where $\mathcal{X} \subseteq \Re^{n}$ is the state space,
$\mathcal{U} \subseteq \Re^{m}$ is the control space,
and
$Q : \mathscr{B}(\mathcal{X}) \times{} \mathcal{X} \times{} \mathcal{U} \rightarrow [0, 1]$
is a stochastic kernel,
which is a Borel-measurable function
that maps a probability measure $Q(\,\cdot\, | \,x, u)$
to each $x \in \mathcal{X}$ and $u \in \mathcal{U}$
on the Borel space $(\mathcal{X}, \mathscr{B}(\mathcal{X}))$.
Further, let $\mathcal{X}$ and $\mathcal{U}$ be compact Borel spaces.
The system evolves over a finite horizon $k \in [0, N]$
with inputs chosen according to a Markov policy
\cite{puterman2014markov, bertsekas2004stochastic}, a sequence
$\pi = \lbrace \pi_{0}, \pi_{1}, \ldots, \pi_{N-1}\rbrace$
of universally-measurable maps
$\pi_{k} : \mathcal{X} \rightarrow \mathcal{U}$.
The set of all Markov control policies $\pi$ is denoted as $\mathcal{M}$.

We define $\mathcal{K}, \mathcal{T} \in \mathscr{B}(\mathcal{X})$
as the \emph{safe set} and \emph{target set}, respectively.
We define the
\emph{terminal-hitting time safety probability}
$r_{x_{0}}^{\pi}(\mathcal{K}, \mathcal{T})$
\cite{summers2010verification}
as the probability that a system $\mathcal{H}$
controlled by a policy $\pi \in \mathcal{M}$
will {reach} $\mathcal{T}$ at $k=N$
while {avoiding} $\mathcal{X}\backslash\mathcal{K}$
for all $k \in [0, N-1]$,
given an initial condition $x_{0} \in \mathcal{X}$.
\begin{align}
  r_{x_{0}}^{\pi}(\mathcal{K}, \mathcal{T}) \triangleq
    \Pr_{x_{0}}^{\pi}
    \lbrace
        \boldsymbol{x}_{N} \in \mathcal{T} \wedge{}
        \boldsymbol{x}_{i} \in \mathcal{K}, \forall i \in [0, N-1]
    \rbrace
  \label{eqn: terminal-hitting probability}
\end{align}
Let $V{}_{k}^{\pi} : \mathcal{X} \rightarrow [0, 1]$
be defined via backward recursion as
\begin{align}
  V{}_{N}^{\pi}(x) &=
	\boldsymbol{1}_{\mathcal{T}}(x)
	\label{eqn: terminal-hitting value N} \\
  V{}_{k}^{\pi}(x) &=
	\boldsymbol{1}_{\mathcal{K}}(x)
  \mathbb{E}_{\boldsymbol{y} \sim Q(\,\cdot\,|\,x, \pi_{k}(x))}
  \left[ V{}_{k+1}^{\pi}(\boldsymbol{y}) \right]

  \label{eqn: terminal-hitting value k}
\end{align}
then $V{}_{0}^{\pi}(x) = r_{x_{0}}^{\pi}(\mathcal{K}, \mathcal{T})$
for every $x_{0} \in \mathcal{X}$.

From \cite[Definition~10]{summers2010verification},
a policy $\pi^{*} \in \mathcal{M}$
is the \emph{maximal reach-avoid policy in the terminal sense}
if and only if
\begin{align}
  r_{x_{0}}^{\pi^{*}}(\mathcal{K}, \mathcal{T}) =
  \sup_{\pi \in \mathcal{M}}
  \lbrace
    r_{x_{0}}^{\pi}(\mathcal{K}, \mathcal{T})
  \rbrace
\end{align}
Let $V{}_{k}^{*} : \mathcal{X} \rightarrow [0, 1]$, $k \in [0, N-1]$
be defined via backward recursion,
initialized with $V{}_{N}^{*}(x) = \boldsymbol{1}_{\mathcal{T}}(x)$, as
\begin{align}
  V{}_{k}^{*}(x) =
  \sup_{u \in \mathcal{U}}
  \left\lbrace
    \boldsymbol{1}_{\mathcal{K}}(x)
    \mathbb{E}_{\boldsymbol{y} \sim Q(\,\cdot\,|\,x,u)}
    \left[ V{}_{k+1}^{*}(\boldsymbol{y}) \right]

  \right\rbrace
  \label{eqn: terminal-hitting optimal value k}
\end{align}
Then, $V{}_{0}^{*}(x) = r_{x_{0}}^{\pi^{*}}(\mathcal{K}, \mathcal{T})$.
If $\pi^{*}$
is such that
\begin{align}
  \pi_{k}^{*}(x) =
  \mathrm{arg}\sup_{u \in \mathcal{U}}
  \left\lbrace
    \boldsymbol{1}_{\mathcal{K}}(x)
    \mathbb{E}_{\boldsymbol{y} \sim Q(\,\cdot \,|\, x, u)}
    \left[ V{}_{k+1}^{*}(\boldsymbol{y}) \right]

  \right\rbrace
  \label{eqn: terminal-hitting optimal control map}
\end{align}
then $\pi^{*}$
is maximal in the terminal sense \cite[Theorem~11]{summers2010verification}.

\subsection{Problem Statement}

Consider a set $\mathcal{S}$ of $M$ samples
of the form $\mathcal{S} = \lbrace (\bar{y}_{i}, \bar{x}_{i}, \bar{u}_{i}) \rbrace_{i = 1}^M$ such that $y_{i}$
is drawn i.i.d. from $Q$
according to $\bar{y}_{i} \sim Q(\,\cdot \,|\, \bar{x}_{i}, \bar{u}_{i})$,
and $\bar{u}_{i} = \pi(\bar{x}_{i})$.
We denote sample vectors with a bar to differentiate them from time-indexed vectors.

\begin{problem}
Without direct knowledge of $Q$,
use samples $\mathcal{S}$
to construct a kernel-based approximation
of \eqref{eqn: terminal-hitting value k}
that converges in probability.
\end{problem}

\begin{problem}
Without direct knowledge of $Q$,
use samples $\mathcal{S}$
to construct a kernel-based approximation
of \eqref{eqn: terminal-hitting optimal value k}
that converges in probability,
in order to compute an approximation of the optimal policy
$\pi^{*}$.
\end{problem}

Using samples $\mathcal{S}$, we employ an approach similar to that in \cite{grunewalder2012modelling}, but that is specific to  the
dynamic program associated with the stochastic reachability problem
for high-dimensional, non-Gaussian systems.
The unique computational efficiencies
afforded by
reproducing kernel Hilbert spaces
transforms computation of
\eqref{eqn: terminal-hitting value k} and
\eqref{eqn: terminal-hitting optimal value k}
into simple matrix operations and inner products.

\section{Kernel Distribution Embeddings for Stochastic~Reachabilty}
\label{section: kernel mean embeddings}

For some set $\mathcal{X}$,
let $\mathscr{H}_{\mathcal{X}}$ denote the unique
reproducing kernel Hilbert space
\cite{scholkopf2001learning}
with the positive definite \cite[Definition~4.15]{steinwart2008support} kernel
$K_{\mathcal{X}} : \mathcal{X} \times{} \mathcal{X} \rightarrow \Re$,
which is a Hilbert space of real-valued functions on $\mathcal{X}$
with inner product
$\langle \,\cdot, \cdot\, \rangle_{\mathscr{H}_{\mathcal{X}}}$
and the induced norm
$\lVert x \rVert_{\mathscr{H}_{\mathcal{X}}} = \left(\langle x, x \rangle_{\mathscr{H}_{\mathcal{X}}} \right)^{1/2}$.
A reproducing kernel Hilbert space has two important properties
\cite{aronszajn1950theory}:
\begin{enumerate}
	\item For any $x, x' \in \mathcal{X}$,
				$K_{\mathcal{X}}(x, \cdot\,) : x' \rightarrow K_{\mathcal{X}}(x, x')$
				is an element of $\mathscr{H}_{\mathcal{X}}$.
	\item An element $K_{\mathcal{X}}(x, x')$ of $\mathscr{H}_{\mathcal{X}}$
				satisfies the \emph{reproducing property}
				such that $\forall f \in \mathscr{H}_{\mathcal{X}}$
				and $x \in \mathcal{X}$,
				\begin{align}
					f(x) &=
					\langle
						K_{\mathcal{X}}(x, \cdot), f(\cdot)
					\rangle_{\mathscr{H}_{\mathcal{X}}}
					\label{eqn: reproducing property} \\

					K_{\mathcal{X}}(x, x') &=
	\langle
		K_{\mathcal{X}}(x, \cdot),
		K_{\mathcal{X}}(x', \cdot)
	\rangle_{\mathscr{H}_{\mathcal{X}}}
				\end{align}

\end{enumerate}

This means that
the evaluation of a function $f \in \mathscr{H}_{\mathcal{X}}$
can be viewed as an inner product
in $\mathscr{H}_{\mathcal{X}}$.
Alternatively, an element $K_{\mathcal{X}}(x, \cdot)$
can be viewed as a nonlinear feature map
$\phi : \mathcal{X} \rightarrow \mathscr{H}_{\mathcal{X}}$,
such that
\begin{align}
	K_{\mathcal{X}}(x, x') =
	\langle
		\phi(x),
		\phi(x')
	\rangle_{\mathscr{H}_{\mathcal{X}}}
\end{align}
Because constructing the feature map $\phi(\,\cdot\,)$
and computing
$\langle \phi(x), \phi(x') \rangle_{\mathscr{H}_{\mathcal{X}}}$
explicitly
can be computationally expensive,
the inner product can be computed using $K_{\mathcal{X}}(x, x')$ directly
for a $K_{\mathcal{X}}$ that is positive definite.
This is known as the \emph{kernel trick} \cite{shawe2004kernel}.

By choosing $K_{\mathcal{X}}$, we effectively choose a basis
to represent the functions in $\mathscr{H}_{\mathcal{X}}$.  With the reproducing property, we can then write
the function $f$ as $f(x) = w^{\top} \phi(x)$, a weighted sum of basis functions for some
possibly infinite-dimensional weight vector $w$.   We wish to solve for the particular $w$ which,
based on the samples $\mathcal S$, minimizes the difference between the observations and the kernel-based estimate.

Let $\mathscr{P}$ denote the set of all probability measures on $\mathcal{X}$.
The \emph{kernel distribution embedding}
\cite{berlinet2011reproducing, smola2007hilbert}
of a probability measure
$\Pr_{\boldsymbol{y}\,|\,x, u} \in \mathscr{P}$,
given by
$\mu_{(x, u)} : \mathscr{P} \rightarrow \mathscr{H}_{\mathcal{X}}$,
is defined as
\begin{align}
	\mu_{(x, u)} \triangleq
	\int_{\mathcal{X}}
	K_{\mathcal{X}}(y, \cdot)

	\Pr_{\boldsymbol{y}\,|\,x, u}
	\lbrace \mathrm{d} y \,|\, x, u \rbrace
\end{align}

Let $\mathscr{H}_{\mathcal{X}}$ denote
the unique RKHS
for the state space $\mathcal{X}$
with the positive definite kernel
$K_{\mathcal{X}} : \mathcal{X} \times{} \mathcal{X} \rightarrow \Re$.
Similarly, let $\mathscr{H}_{\mathcal{X} \times \mathcal{U}}$ denote
the RKHS
for $\mathcal{X} \times \mathcal{U}$
with the positive definite kernel
$K_{\mathcal{X} \times \mathcal{U}} : (\mathcal{X}, \mathcal{U}) \times{} (\mathcal{X}, \mathcal{U}) \rightarrow \Re$.
We define the \emph{conditional distribution embedding} of the stochastic kernel
$Q \in \mathscr{P}$ as $\mu_{(x, u)}$.
Then, the expectation of $f$ with respect to the probability measure $Q$ is given by
\begin{align}
	\langle \mu_{(x,u)}, f \rangle_{\mathscr{H}_{\mathcal{X}}} =
	\mathbb{E}_{\boldsymbol{y} \sim Q(\,\cdot\,|\,x, u)}
	[ f(\boldsymbol{y}) ]
	\label{eqn: conditional distribution embedding inner product}
\end{align}
This means we can evaluate the expectation of a function
with respect to $Q$ as an inner product
in $\mathscr{H}_{\mathcal{X}}$.

We can construct an estimate $\bar{\mu}_{(x, u)}$ of $\mu_{(x, u)}$
\cite{smola2007hilbert}
from samples $\mathcal{S}$
to approximate
\eqref{eqn: conditional distribution embedding inner product},
\begin{align}
	\langle \bar{\mu}_{(x, u)}, f \rangle_{\mathscr{H}_{\mathcal{X}}} \approx
	\mathbb{E}_{\boldsymbol{y} \sim Q(\,\cdot \,|\, x, u)}
	\left[
		f(\boldsymbol{y})
	\right]
	\label{eqn: conditional distribution embedding inner product estimate}
\end{align}
According to the Riesz representation theorem
\cite{micchelli2005learning},
the element $\bar{\mu}_{(x, u)}$
can be viewed as the solution to a regularized least-squares problem
that minimizes the error of the expectation operator over the samples \cite{micchelli2005learning, lever2012conditional},
\begin{align}

    \min_{\bar{\mu}} \left\lbrace
		\frac{1}{M}
    \sum_{i=1}^{M} \lVert K_{\mathcal{X}}(\bar{y}_{i}, \cdot) - \bar{\mu}_{(\bar{x}_{i}, \bar{u}_{i})} \rVert_{\mathscr{H}_{\mathcal{X}}}^{2} + \lambda \lVert \bar{\mu} \rVert_{\mathscr{H}_{\Upsilon}}^{2} \right\rbrace
		\label{eqn: regularized least-squares}
\end{align}
where $\mathscr{H}_{\Upsilon}$ is a vector-valued RKHS \cite{micchelli2005learning}.
The solution is unique and has the form
\begin{align}
	\bar{\mu}_{(x, u)} =
	\eta
	\Phi^{\top}
	(G + \lambda M I)^{-1}
	\Psi
	K_{\mathcal{X} \times \mathcal{U}}((x, u), \cdot)

	\label{eqn: conditional distribution embedding estimate}
\end{align}
where $\lambda$ is a regularization parameter to avoid overfitting,
$\eta$ is a normalizing constant,
and $\Phi$, $\Psi$, and $G$ are defined as
\begin{align}
	\Phi &= [
		K_{\mathcal{X}}(\bar{y}_{1}, \cdot),
		\ldots,
		K_{\mathcal{X}}(\bar{y}_{M}, \cdot)
	]^{\top}
	\label{eqn: phi vector} \\
	\Psi &= [
		K_{\mathcal{X} \times \mathcal{U}}((\bar{x}_{1}, \bar{u}_{1}), \cdot),
		\ldots,
		K_{\mathcal{X} \times \mathcal{U}}((\bar{x}_{M}, \bar{u}_{M}), \cdot)
	]^{\top}
	\label{eqn: psi vector} \\
	G &= \Psi \Psi^{\top}
	\label{eqn: gram matrix}
\end{align}
By the reproducing property
of $K_{\mathcal{X}}$ in $\mathscr{H}_{\mathcal{X}}$,
$\forall f \in \mathscr{H}_{\mathcal{X}}$,
we can rewrite
\eqref{eqn: conditional distribution embedding inner product estimate} as
\begin{align}
  \langle \bar{\mu}_{(x, u)}, f \rangle_{\mathscr{H}_{\mathcal{X}}} &=
		\boldsymbol{f}^{\top} \beta(x, u)
	\label{eqn: conditional distribution embedding inner product estimate as sum}
\end{align}
where $\boldsymbol{f} = [f(\bar{y}_{1}), \ldots, f(\bar{y}_{M})]^{\top}$,
and $\beta(x, u) \in \Re^{M}$ is a vector of coefficients,
\begin{align}

	\beta(x, u) &=
	\eta
	(G + \lambda M I)^{-1}
	\Psi
	K_{\mathcal{X} \times \mathcal{U}}((x, u), \cdot)
	\label{eqn: compute beta}
\end{align}
\emph{
This means an approximation of the value function expectation
$\mathbb{E}_{\boldsymbol{y} \sim Q(\,\cdot \,|\, x, \pi_{k}(x))} \left[ V{}_{k+1}^{\pi}(\boldsymbol{y}) \right]$
in \eqref{eqn: terminal-hitting value k}
can be evaluated as a linear operation
in $\mathscr{H}_{\mathcal{X}}$.
}

\subsection{Terminal-Hitting Time Problem}
\label{section: terminal-hitting time problem}

With the conditional distribution embedding $\mu_{(x, u)}$,
the value functions in \eqref{eqn: terminal-hitting value k}
can be written as
\begin{align}
	V{}_{k}^{\pi}(x) =
	\boldsymbol{1}_{\mathcal{K}}(x)
	\langle
		\mu_{(x, \pi_{k}(x))}, V{}_{k+1}^{\pi}
	\rangle_{\mathscr{H}_{\mathcal{X}}}
	\label{eqn: terminal-hitting value k as inner product}
\end{align}
With the estimate $\bar{\mu}_{(x, u)}$
\eqref{eqn: conditional distribution embedding estimate},
we define the approximate value functions
$\overbar{V}_{k}^{\pi} : \mathcal{X} \rightarrow [0, 1]$, $k \in [0, N-1]$, as
\begin{align}
	\overbar{V}_{k}^{\pi}(x) =
	\boldsymbol{1}_{\mathcal{K}}(x)
	\langle
		\bar{\mu}_{(x, \pi_{k}(x))}, V{}_{k+1}^{\pi}
	\rangle_{\mathscr{H}_{\mathcal{X}}}
	\label{eqn: terminal-hitting approximate value k as inner product}
\end{align}
such that $\overbar{V}_{k}^{\pi}(x) \approx	V{}_{k}^{\pi}(x)$.
Thus, an approximation for the reach-avoid probability
$r_{x_{0}}^{\pi}(\mathcal{K}, \mathcal{T})$
can be computed via the backward recursion
described in \algorithmref{algo: backward recursion},
such that
\begin{align}
	\overbar{V}_{0}^{\pi}(x)
	\approx r_{x_{0}}^{\pi}(\mathcal{K}, \mathcal{T})
	\label{eqn: terminal-hitting approximate value 0}
\end{align}

\begin{algorithm}[!t]
	\caption{Value Function Estimate}
	\label{algo: backward recursion}
	\textbf{Input}:
	samples $\mathcal{S}$ drawn i.i.d. from $Q$,

	policy $\pi$,
	horizon $N$
	\\
	\textbf{Output}:
	value function estimate $\overbar{V}_{0}^{\pi}(x)$

	\begin{algorithmic}[1]
		\State Compute $\Phi$ and $\Psi$ using $\mathcal{S}$ from
		\eqref{eqn: phi vector} and
		\eqref{eqn: psi vector}
		\State $G \gets \Psi \Psi^{\top}$
		\State Initialize $\overbar{V}_{N}^{\pi}(x) \gets \boldsymbol{1}_{\mathcal{T}}(x)$
		\For{$k \gets N-1$ to $0$}

			\State Compute $\beta(x, \pi_{k}(x))$ from \eqref{eqn: compute beta}
			\State $

            \mathcal{Y}
			\gets
			[\overbar{V}_{k+1}^{\pi}(\bar{y}_{1}), \ldots, \overbar{V}_{k+1}^{\pi}(\bar{y}_{M})]^{\top}$
			\State
				$
					\overbar{V}_{k}^{\pi}(x) \gets
						\boldsymbol{1}_{\mathcal{K}}(x)

				\mathcal{Y}^{\top}
							\beta(x, \pi_{k}(x))
				$

		\EndFor
		\State Return $\overbar{V}_{0}^{\pi}(x) \approx r_{x_{0}}^{\pi}(\mathcal{K}, \mathcal{T})$
	\end{algorithmic}

\end{algorithm}

We now seek to characterize
the quality of the approximation
and the conditions for its convergence.
As in \cite{sriperumbudur2010hilbert, gretton2007kernel, gretton2012kernel},
we define a pseudometric
that characterizes the accuracy of the estimate $\bar{\mu}_{(x,u)}$.
\begin{definition}[Distance Pseudometric]
The distance pseudometric in $\mathscr{H}_{\mathcal{X}}$ between
the conditional distribution embedding
$\mu_{(x, u)} \in \mathscr{H}_{\mathcal{X}}$
and the estimate
$\bar{\mu}_{(x, u)} \in \mathscr{H}_{\mathcal{X}}$ is defined as
$\lVert \mu_{(x, u)} - \bar{\mu}_{(x, u)} \rVert_{\mathscr{H}_{\mathcal{X}}}$.
\end{definition}

It is shown in \cite{fukumizu2008kernel}
that if $K_{\mathcal{X}}$ is a \emph{characteristic}, bounded kernel,
then
$\lVert \mu_{(x, u)} - \bar{\mu}_{(x, u)} \rVert_{\mathscr{H}_{\mathcal{X}}} = 0$
if and only if $\mu_{(x, u)} = \bar{\mu}_{(x, u)}$.
A kernel is characteristic if the kernel embedding is injective,
meaning the embeddings for any two different conditional distributions
are represented by different elements
within the RKHS.
Thus, as
$\lVert \mu_{(x, u)} - \bar{\mu}_{(x, u)} \rVert_{\mathscr{H}_{\mathcal{X}}}$
converges
\cite{grunewalder2012modelling, smola2007hilbert},
the estimate converges in probability
to the conditional distribution embedding
within $\mathscr{H}_{\mathcal{X}}$.

\begin{lemma}\cite[Lemma~2.2]{grunewalder2012modelling}
	\label{lemma: conditional distribution embedding convergence}

	For any $\varepsilon > 0$,
	if the regularization parameter $\lambda$ in
	\eqref{eqn: regularized least-squares}
	is chosen such that
	$\lambda \rightarrow 0$
	and $\lambda^{3} M \rightarrow \infty$, and if $\mathcal{X}$ is bounded and $K_{\mathcal{X}}$ is strictly positive definite, then
	\begin{align}
		\Pr_{\mathcal{S} \sim Q}
		\left\lbrace
			\sup_{(x, u) \in \mathcal{X} \times \mathcal{U}}
			\lVert
				\mu_{(x, u)} - \bar{\mu}_{(x, u)}
			\rVert_{\mathscr{H}_{\mathcal{X}}}
			> \varepsilon
		\right\rbrace \rightarrow 0
		\label{eqn: conditional distribution embedding convergence}
	\end{align}

\end{lemma}

\begin{proposition}[Value Function Convergence]
	\label{prop: value function expectation convergence}

	For any $\varepsilon > 0$,
	if the regularization prameter $\lambda$ in
	\eqref{eqn: regularized least-squares}
	is chosen such that
	$\lambda \rightarrow 0$
	and $\lambda^{3} M \rightarrow \infty$, and if $\mathcal{X}$ is bounded and $K_{\mathcal{X}}$ is strictly positive definite,
	$\vert V{}_{k}^{\pi}(x) - \overbar{V}_{k}^{\pi}(x) \vert$
	converges in probability.

\end{proposition}

\begin{proof}
	By subtracting
	\eqref{eqn: terminal-hitting approximate value k as inner product} from
	\eqref{eqn: terminal-hitting value k as inner product}, and using the parallelogram law, we define
	the absolute value function error $\mathcal{E}_{k}(x)$ at time $k$,
	\begin{align}
		\mathcal{E}_{k}(x) &\triangleq
		\vert
			V{}_{k}^{\pi}(x) - \overbar{V}_{k}^{\pi}(x)
		\vert \\

&= \boldsymbol{1}_{\mathcal{K}}(x)
		\vert
			\langle
				\mu_{(x,\pi_{k}(x))} - \bar{\mu}_{(x,\pi_{k}(x))}, V{}_{k+1}^{\pi}
			\rangle_{\mathscr{H}_{\mathcal{X}}}
		\vert
		\label{eqn: value function absolute difference}
	\end{align}
	We can rewrite
	\eqref{eqn: value function absolute difference}
	using Cauchy-Schwarz to obtain
	\begin{align}
		\mathcal{E}_{k}(x)

		&\leq \boldsymbol{1}_{\mathcal{K}}(x)
		\lVert
			V{}_{k+1}^{\pi}
		\rVert_{\mathscr{H}_{\mathcal{X}}}
		\lVert
			\mu_{(x,\pi_{k}(x))} - \bar{\mu}_{(x,\pi_{k}(x))}
		\rVert_{\mathscr{H}_{\mathcal{X}}}
	\end{align}
	Since
	$\lVert \mu_{(x,\pi_{k}(x))} - \bar{\mu}_{(x,\pi_{k}(x))} \rVert_{\mathscr{H}_{\mathcal{X}}}$
	converges in probability
	according to
	\lemmaref{lemma: conditional distribution embedding convergence},
	$\vert V{}_{k}^{\pi}(x) - \overbar{V}_{k}^{\pi}(x) \vert$
	also converges in probability
	with the probabilistic error bound $\varepsilon$.
\end{proof}
Using this, the value function approximation in
\eqref{eqn: terminal-hitting approximate value k as inner product}
converges in probability
for some probabilistic error bound $\varepsilon$
as the number of samples increases.

\begin{corollary}
	\label{cor: value function expectation convergence cor1}
	For any $\varepsilon > 0$,
	the error in the reach-avoid probability
	computed using \algorithmref{algo: backward recursion}

	converges in probability to
	\begin{align}
		\vert
			V{}_{0}^{\pi}(x) - \overbar{V}_{0}^{\pi}(x)
		\vert \leq N\varepsilon
	\end{align}
\end{corollary}
\begin{proof}
	By subtracting
	\eqref{eqn: terminal-hitting approximate value k as inner product} from
	\eqref{eqn: terminal-hitting value k as inner product}, we obtain
	the absolute value function error $\mathcal{E}_{N-1}(x)$ at time $k = N-1$,
	\begin{align}
		\mathcal{E}_{N-1}(x) =
		\vert
			V{}_{N-1}^{\pi}(x) - \overbar{V}_{N-1}^{\pi}(x)
		\vert

		\label{eqn: appendix cor1 value function k error}
	\end{align}
	Using \propositionref{prop: value function expectation convergence},
    if the error in the approximate value function converges in probability to at most $\varepsilon$,
	then the error in \eqref{eqn: appendix cor1 value function k error} converges in probability to $\varepsilon$,
	i.e. $\mathcal{E}_{N-1}(x) \leq \varepsilon$.

	Because the error
	in the approximate value function
	for $k = N-1$
	converges in probability to $\varepsilon$,
	then by approximating and recursively substituting
	$\overbar{V}_{k}^{\pi}(x)$ for $k < N-1$,
	the error at time $k$ converges in probability to $(N-k)\varepsilon$.

	Thus, by induction, the error
	obtained by the backward recursion in
	\algorithmref{algo: backward recursion}
	converges in probability to at most $N\varepsilon$,
	\begin{align}

		\vert
			V{}_{0}^{\pi}(x) - \overbar{V}_{0}^{\pi}(x)
		\vert \leq N\varepsilon
	\end{align}
	which concludes the proof.
\end{proof}

\subsection{Maximal Reach-Avoid Policy in the Terminal Sense}

As in \eqref{eqn: terminal-hitting value k as inner product},
we write the optimal value functions $V{}_{k}^{*}$
from \eqref{eqn: terminal-hitting optimal value k}
using the conditional distribution embedding $\mu_{(x, u)}$.
\begin{align}
	V{}_{k}^{*}(x) =
	\sup_{u \in \mathcal{U}}
	\left\lbrace
		\boldsymbol{1}_{\mathcal{K}}(x)
		\langle
			\mu_{(x, u)}, V{}_{k+1}^{*}
		\rangle_{\mathscr{H}_{\mathcal{X}}}
	\right\rbrace
	\label{eqn: terminal-hitting optimal value k as inner product}
\end{align}
With the estimate $\bar{\mu}_{(x, u)}$
from \eqref{eqn: conditional distribution embedding estimate},
we define the approximate optimal value functions
$\overbar{V}_{k}^{*} : \mathcal{X} \rightarrow [0, 1]$, $k \in [0, N-1]$
\begin{align}
	\overbar{V}_{k}^{*}(x) =
	\sup_{u \in \mathcal{U}}
	\left\lbrace
	\boldsymbol{1}_{\mathcal{K}}(x)
	\langle
		\bar{\mu}_{(x, u)}, V{}_{k+1}^{*}
	\rangle_{\mathscr{H}_{\mathcal{X}}}

	\right\rbrace
	\label{eqn: terminal-hitting approximate optimal value k as inner product}
\end{align}
such that $\overbar{V}_{k}^{*}(x) \approx V{}_{k}^{*}(x)$.
If
$\bar{\pi}_{k}^{*} : \mathcal{X} \rightarrow \mathcal{U}$
is such that
\begin{align}
	\bar{\pi}_{k}^{*}(x) =
	\mathrm{arg}\sup_{u \in \mathcal{U}}
	\left\lbrace
		\boldsymbol{1}_{\mathcal{K}}(x)
		\langle
			\bar{\mu}_{(x, u)}, V{}_{k+1}^{*}
		\rangle_{\mathscr{H}_{\mathcal{X}}}
	\right\rbrace
	\label{eqn: approximate maximal reach-avoid policy in the terminal sense}
\end{align}
then
$\bar{\pi}^{*} = \lbrace \bar{\pi}_{0}^{*}, \bar{\pi}_{1}^{*}, \ldots \rbrace$
is the approximate maximal reach-avoid policy in the terminal sense.
The approximate optimal reach-avoid probability
under policy $\bar{\pi}^{*}$
initialized with $\overbar{V}_{k}^{*}(x) = \boldsymbol{1}_{\mathcal{T}}(x)$,
is described in
\algorithmref{algo: backward recursion} as
\begin{align}
  r_{x_{0}}^{*}(\mathcal{K}, \mathcal{T}) \approx
	\overbar{V}_{0}^{*}(x).
	\label{eqn: terminal-hitting approximate optimal value 0}
\end{align}

\section{Numerical Results}
\label{section: numerical results}

\begin{table*}[ht!]
	\caption{Computation Time and Parameters}
	\label{table: computation time}
	\centering
	\normalsize
	\renewcommand{\arraystretch}{1.1}
	\begin{tabular*}{\textwidth}{
		@{\extracolsep{\fill}} l
		@{\hspace{0.5\tabcolsep}} c
		@{\hspace{0.5\tabcolsep}} c
		@{\hspace{0.5\tabcolsep}} c
		@{\hspace{0.5\tabcolsep}} c
		@{\hspace{0.5\tabcolsep}} c
		}
		\toprule
		& Number of
		& Number of
		&
		& Dynamic
		& Chance-Constrained
		\\
		  System
		& Samples [$M$]
		& Evaluation Points
		& \algorithmref{algo: backward recursion}
		& Programming
		& Open \cite{vinod2019affine}
		\\
		\midrule
		Double Integrator
		& $1024$
		& $10201$
		& $0.43$ s
		& $31.76$ s
		& $24.36$ s
		\\
		CWH
		& $883$
		& $10201$
		& $0.48$ s
		& --
		& $34.51$ s
		\\
		\midrule
		10000-D Integrator
		& $1024$
		& $1$
		& $30.53$ s
		& --
		& --
		\\
		\bottomrule
	\end{tabular*}
\end{table*}

We implemented \algorithmref{algo: backward recursion} on two
well-known stochastic systems for the purpose of validation and error analysis.
We first generated $M = 1024$ samples via simulation,
then presumed no knowledge of the system dynamics
or the stochastic disturbance in computing
$r_{x_{0}}^{\pi}(\mathcal{K}, \mathcal{T})$ via Algorithm \ref{algo: backward recursion}.
We chose a
Gaussian radial basis function kernel,
 $K(x, x') = \exp ( - \lVert x - x' \rVert_{2}^{2} / 2\sigma^2 )$,
with $\sigma = 0.1$,
and regularization parameter $\lambda = 1$.
All computations were done in Matlab on a
2.7GHz Intel Core i7 CPU with 16 GB RAM.
Code to generate all figures is available at
\textit{https://github.com/unm-hscl/ajthor-LCSS-2019}.

\subsection{$n$-D Stochastic Chain of Integrators}

We consider a $n$-D stochastic chain of integrators
\cite{vinod2017scalable}, in which the input appears at the $n^{\mathrm{th}}$ derivative
and each element of the state vector is the discretized integral of the element that follows it. The dynamics with sampling time $T$ are given by:
\begin{align}
    \boldsymbol{x}_{k+1} =
    \begin{bmatrix}
      1 & T & \cdots & \frac{T^{n-1}}{(n-1)!} \\
      0 & 1 & \cdots & \frac{T^{n-2}}{(n-2)!} \\
      \vdots & \vdots & \ddots & \vdots \\
      0 & 0 & 0 & 1
    \end{bmatrix}
    \boldsymbol{x}_{k} +
    \begin{bmatrix}
      \frac{T^{n}}{n!} \\
      \frac{T^{n-1}}{(n-1)!} \\
      \vdots \\
      T
    \end{bmatrix}
    u_{k} +
    \boldsymbol{w}_{k}
  \end{align}
with i.i.d. disturbance $\boldsymbol{w}_{k}$
defined on the probability space
$(\mathcal{W}, \mathscr{B}(\mathcal{W}), \Pr_{\boldsymbol{w}})$.
We consider two distributions:
1) a Gaussian distribution with variance $\Sigma = 0.01I$
such that $\boldsymbol{w}_{k} \sim \mathcal{N}(0, \Sigma)$, and
2) a beta distribution such that $\boldsymbol{w}_{k} \sim \textnormal{Beta}(\alpha, \beta)$, with PDF
	$f(x \,|\, \alpha, \beta) =
	\frac{
		\Gamma(\alpha + \beta)
	}{
		\Gamma(\alpha) \Gamma(\beta)
	}
	x^{\alpha-1} (1-x)^{\beta-1}$,
described in terms of Gamma function $\Gamma$ and
positive shape parameters $\alpha = 0.5$, $\beta = 0.5$.
The control policy $\pi$ is
$\pi_{0}(x) = \pi_{1}(x) = \ldots = 0$.
The target set and safe set are $\mathcal{T}, \mathcal{K} = [-1, 1]^{n}$.

For a system with $n=2$, we compute
approximate safety probabilities
using \algorithmref{algo: backward recursion}
for time horizon $N=3$ with the Gaussian distribution in Fig. \ref{fig: main_figure_1}(a) and
the beta
distribution in Fig. \ref{fig: main_figure_2}(a).
We compared the RKHS solution for the Gaussian distribution with a dynamic programming solution
implemented in \cite{vinod2019sreachtools}.
The absolute error (\ref{eqn: value function absolute difference})
(Fig. \ref{fig: main_figure_1}(b)) has a maximum value of 0.074.
We consider the region strictly within $\mathcal{K}$ to account for Matlab rounding errors.
The error is highest along the ridges on the upper right and lower left corners.
Fig. \ref{fig: dimensionality_time} shows that as the number of samples increases, the error decreases, as expected.

For a 10,000-dimensional system, which is beyond the
computational capabilities of any existing methods for stochastic reachability,
computation of \eqref{eqn: terminal-hitting approximate value 0} took 30.53 seconds for $x_0 = 0$ (Table \ref{table: computation time}).
Computation time, evaluated from the
the same initial condition for all systems of dimension 2 through 10,000 (Fig. \ref{fig: dimensionality_time}),
appears to increase linearly
because computation of the norm in the Gaussian kernel function
scales linearly with state dimension.
Note that the structure of the system dynamics plays no role in the computational complexity,
as the structure of $G$ in Algorithm~\ref{algo: backward recursion} does not depend on the dynamics.

\begin{figure}[t!]
    \centering
    \includegraphics{main_figure_1}
    \caption{
    (a) Approximate safety probabilities computed using \algorithmref{algo: backward recursion} for a double integrator at $k=0$ for $N=3$.
	(b) Absolute error $\vert V_{0}^{\pi}(x) - \overbar{V}_{0}^{\pi}(x) \vert$
	between the dynamic programming solution
	and \algorithmref{algo: backward recursion}
	at $k=0$ for $N=3$.}
    \label{fig: main_figure_1}
\end{figure}

\begin{figure}[!t]
    \centering
	\includegraphics{dim_time}
	\caption{
		(top) System dimensionality [$n$] vs. average computation time [$s$] for an $n$-D stochastic integrator system.
		(bottom) Number of samples [$M$] vs. maximum absolute error $\vert V_{0}^{\pi}(x) - \overbar{V}_{0}^{\pi}(x) \vert$ for time steps $k = 2$ to $k = 0$ for $N = 3$.
	}
	\label{fig: dimensionality_time}
\end{figure}

\subsection{Clohessy-Wiltshire-Hill System}

Lastly, we considered the more realistic example of
spacecraft rendezvous and docking, in which
a spacecraft must dock with another spacecraft
while remaining within a line of sight cone.
The Clohessey-Wiltshire-Hill dynamics,
\begin{align}
		\ddot{x} - 3\omega x - 2 \omega \dot{y} &= F_{x}/m_{d} &
		\ddot{y} + 2 \omega \dot{x} &= F_{y}/m_{d}
\end{align}
with state $\boldsymbol{z} = [x, y, \dot{x}, \dot{y}] \in \mathcal{X} \subseteq \Re^{4}$,
 input
 $u = [F_{x}, F_{y}] \in \mathcal{U} \subseteq \Re^{2}$,
where $\mathcal{U} = [-0.1, 0.1] \times [-0.1, 0.1]$, and parameters $\omega$, $m_d$
can be written as a discrete-time LTI system
$\boldsymbol{z}_{k+1} = A \boldsymbol{z}_{k} + B u_{k} + \boldsymbol{w}_{k}$
with an additive Gaussian disturbance \cite{lesser2013stochastic}
with variance $\Sigma = \textnormal{diag}(1 \times 10^{-4}, 1 \times 10^{-4}, 5 \times 10^{-8}, 5 \times 10^{-8})$
such that $\boldsymbol{w}_{k} \sim \mathcal{N}(0, \Sigma)$.
The target set and safe set are defined
as in \cite{lesser2013stochastic}:
\begin{align}
	\mathcal{T} &=
	\lbrace
		\boldsymbol{z} \in \Re^{4} :
		\begin{aligned}[t]
				& \vert \boldsymbol{z}_{1} \vert \leq 0.1,
					-0.1 < \boldsymbol{z}_{2} < 0, \\
				& \vert \boldsymbol{z}_{3} \vert \leq 0.01,
					\vert \boldsymbol{z}_{4} \vert \leq 0.01 \rbrace
		\end{aligned} \\
	\mathcal{K} &=
	\lbrace
		\boldsymbol{z} \in \Re^{4} 	: \vert \boldsymbol{z}_{1} \vert < \vert \boldsymbol{z}_{2} \vert,
											\vert \boldsymbol{z}_{3} \vert \leq 0.05,
											\vert \boldsymbol{z}_{4} \vert \leq 0.05
	\rbrace
\end{align}
We generate samples using \cite{vinod2019sreachtools}
with a chance-affine controller.

Fig. \ref{fig: main_figure_2}(b) shows the approximate safety probabilities for time horizon $N = 5$, with a precomputed safety controller from \cite{vinod2019affine, vinod2019sreachtools}.
The safety probabilities for the entire region were computed in 0.48 seconds (Table \ref{table: computation time}), almost two orders of magnitude less than the chance constrained approach \cite{vinod2019affine}, which computes the set of initial conditions
where the safety probability is above a certain threshold (0.8 in this case).

\begin{figure}[t!]
    \centering
    \includegraphics{main_figure_2}
    \caption{(a) Approximate safety probabilities computed using \algorithmref{algo: backward recursion} for a double integrator
	with a
	$\textnormal{Beta}(0.5, 0.5)$
	disturbance at $k=0$ for $N=1$.
	(b) Terminal-hitting safety probabilities for a CWH system
	with a Gaussian disturbance
	and a chance-affine controller
	at $k=0$ for $N=5$.}
    \label{fig: main_figure_2}
\end{figure}

\subsection{Sample Size and Parameter Tuning}

The number of samples used to create the estimate $\bar{\mu}_{(x, u)}$
is the most significant computational bottleneck for
\algorithmref{algo: backward recursion} and is generally $\mathcal{O}(M^3)$.
As the number of samples
increases,
the absolute error in the approximate safety probabilities decreases.
However, methods have been developed recently to
reduce the computational complexity \cite{rahimi2008random, le2013fastfood, lever2012conditional}
to $\mathcal{O}(M \log M)$.
Additionally, for high-dimensional systems,
the number of samples needed to fully characterize the
system dynamics and disturbance
increases as the system dimensionality increases,
which is prohibitive for analysis over large regions of the state space.
However, due to the sample-based nature of
\algorithmref{algo: backward recursion}, we can choose samples
within a local region of interest
in order to approximate the safety probabilities.

The kernel bandwidth parameter and the regularization parameter
are tunable parameters which can affect
the quality of the estimate obtained
using \algorithmref{algo: backward recursion}.
A cross-validation scheme to empirically
choose these parameters is presented in
\cite[Section~6]{micchelli2005learning}.

\section{Conclusions \& Future Work}

We present a sample-based method to compute the stochastic reachability probability measure
for Markov control systems
with arbitrary disturbances that does not require a known model of the transition kernel.
Our approach employs efficient computation associated with a reproducing kernel Hilbert space
to approximate conditional distributions via simple matrix operations.
The method is demonstrated on a 10,000-dimensional integrator as well as a realistic model of
relative spacecraft motion.
We plan to extend this
to sample-based controller synthesis, with application to
systems with autonomous and human elements.

\bibliographystyle{IEEEtran}


\section{Introduction}

The most important experiments which have estimated the amount of Dark
Matter in the present universe are Cosmic Microwave Background Radiation
(CMBR) experiments, discussed in the section 2.

There have been a number of theoretical models for the creation of Dark
Matter, which is reviewed in section 3. It is almost certain that sterile
nuetrinos are part of Dark Matter.

Experiments detecting sterile nuetrinos via neutrino oscillation and
a theoretical study of neutrino oscillations with 3 active and 3 sterile
neutrinos with the present results are discussed in section 4. Also a recent
search for sub-Gev Dark Matter by the MiniBooNE-DM Collaboration is briefly
discussed in section 4.

Neutrinos are fermions with quantum spin 1/2. It is possible that some Dark
Matter particle are vector bosons with quantum spin 1, like the photon. A
possibe Dark Photon Search is discussed in section 5. Future experimental
detection of Dark Matter is reviewed in section 6.

\section{Cosmic Microwave Background Radiation (CMBR)}

There have been many CMBR experiments, such as Refs\cite{cobe96,acbar08,quad10,
wmap13} which have estimated the total density of the present universe,
Dark Matter density, dark energy density, etc.

With $\Omega$ the density of the Universe and $\Omega=1.0$ for a flat Universe,
results from CMBR observations are:
\begin{eqnarray}
\label{CMBR-results}
\Omega &=& 1.0023^{+.0056}_{-.0054} \nonumber \\
{\rm Dark\;Energy\;Density\;(vacuum\;energy)}&\simeq& 0.73 \nonumber \\
{\rm Dark\;Matter\;Density}&\simeq& 0.23 \nonumber \\
{\rm Baryon (Normal Matter)\;Density} &\simeq& 0.04 \nonumber \\
{\rm Age\;of\;the\;Universe}&\simeq& 1.37 {\rm\;\;billion\;years} \; .
\end{eqnarray}
Therefore about 23 \% of the Universe is Dark Matter.

About 73 \% of the universe is dark energy, which is anti-gravity. Dark
Energy (Quintessence) is anti-gravity and produced inflation at a very
early time, which is why we now have an almost homogeneous universe.

Only about 4 \% of the universe is normal matter.
\newpage
\section{Creation of Dark Matter}

There have been many studies of the origin of Dark Matter mass. In order
to explain why the universe is highly homongeneous Guth\cite{guth81} introduced
what is now called the dark energy or quintessence field to produce inflation.
Peebles and Ratra\cite{pr88} studied the quintessence field. Farrar and
Peebles\cite{fp04} intersecting Dark Matter and the Dark Energy (quintessence)
field.

A recent study estimated Dark Matter mass created via Dark Energy
interaction during the Cosmological Electroweak Phase Transition
 (EWPT)\cite{kc14} at a time $t_{EWPT} = 10^{-11}$ s.
Using the Dark Energy Lagrangian ${\cal L}^{DE}$
and Dark Mass-Dark Energy Lagrangian ${\cal L}^{DM-DE}$
\begin{eqnarray}
\label{L}
  {\cal L}^{DE}&=& \frac{1}{2}\partial_\nu \Phi_q \partial^\nu \Phi_q
-V(\Phi_q) \nonumber \\
{\cal L}^{DM-DE}&=& g_D \bar{\psi}^{DM} \Phi_q \psi^{DM} \; ,
\end{eqnarray}
with $\Phi_q$ the quintessence field and $V(\Phi_q)$ obtained from
Refs\cite{pr88,fp04}. Depending on the choice of parameters Ref\cite{kc14}
estimated Dark Matter mass $M_{DM}(EWPT)$ at the present time $t=t_{now}$  as
\begin{eqnarray}
\label{DM-DE-EWPT}
M_{DM}(EWPT) &\simeq& {\rm few\;GeV\;\;to\;140\;GeV} \ .
\end{eqnarray}
With this range of values Dark Matter particles might be detected in the
near future.

More recently a study estimated Dark Matter mass created via Dark Energy
interaction during the Cosmological Quantum Chromodynamics Phase Transition
(QCDPT)\cite{kd19} using the Lagrangians in Eq(\ref{L}) at a time $t_{QCDPT} =
 10^{-4}$ s, with the result at $t=t_{now}$
\begin{eqnarray}
\label{DM-DE-QCDPT}
M_{DM}(QCDPT) &\simeq& {\rm 0.5\;\;to\;3.5\;TeV} \ ,
\end{eqnarray}
which is more than an order of magnitude larger than $M_{DM}(EWPT)$ as
$t_{QCDPT}\gg t_{EWPT}$.

\section{Neutrino Oscillations and Sterile Neutrino Parameters}

Sterile neutrinos are a well-known source of Dark Matter.
For many years the MiniBooNE Collaboration has studied neutrino oscillations.
The MiniBooNE  Collaboration\cite{mini13} using $\bar{\nu}_\mu \rightarrow
\bar{\nu}_e$ oscillations detected a sterile neutrino $\nu_4$
and estimated the mass difference between this sterile neutrino and a standard
neutrino as
\begin{eqnarray}
\label{sterilemass}
\Delta m^2&=& m_4^2-m_1^4 \simeq 0.06 {\rm (eV)}^2 \; ,
\end{eqnarray}
which is too small for sterile neutrino $\nu_4$ to be a WIMP (Weakly Interacting
Massive Particle).

In a calculation with three active and three sterile neutrinos\cite{lsk16},
an extension of work with six neutrinos\cite{tg07},
$\mathcal{P}(\nu_\mu\rightarrow \nu_e)$, the transition probability for a
muon neutrino to oscillate to an electron neutrino, was estimated.

With the notation for the three active neutrinos $\nu_e=\nu_1,\nu_\mu=\nu_2,
\nu_\tau=\nu_3$ and the three sterile neutrinos are $\nu_4,\nu_5,\nu_6$
\begin{eqnarray}
\label{Pue-1}
 \mathcal{P}(\nu_\mu \rightarrow\nu_e) &=& Re[\sum_{i=1}^{6}\sum_{j=1}^{6}
U_{1i}U^*_{1j}U^*_{2i}U_{2j} e^{-i(\delta m_{ij}^2/E)L}] \; ,
\end{eqnarray}
where neutrinos with flavors $\nu_e,\nu_\mu,\nu_\tau$ and three sterile
neutrinos, $\nu_{s_1},\nu_{s_2},\nu_{s_3}$ are related to neutrinos with
definite mass by
\begin{eqnarray}
\label{f-mrelation}
      \nu_f &=& U\nu_m \; ,
\end{eqnarray}
where $U$ is a 6x6 matrix. In Ref\cite{lsk16} the sterile-active neutrino
mixing angles were calculated, which should be useful for future searches
for sterile neutrinos as Dark Matter particles.

Recently the MiniBooNE-DM Collaboration has carried out a search for sub-Gev
Dark Matter at the Fermilab\cite{mbde18}, but no excess from background
predictions was observed.

\section{Dark Photon Search}

There is a worldwide race in the search for Dark Photons.

A Dark Photon would
be a vector boson (quantum spin=1), while a neutrino is fermion (quantum
spin=1/2, so Dark Photons would be a very different kind of Dark Matter
particle than a sterile neutrino.

The Dark Photon Search is described in Ref\cite{p25}. The LANL(Los Alamos
National Laoratory), P-25, PHENIX team is planning to search
for Dark Photons at Fermilab. The Direct Search for Dark Photons at Fermilab
is reviewed in Ref\cite{kl16}, which can be found by clicking on Ref\cite{kl16}
in ``talks and presentations'' in Ref\cite{p25}. Since the detection of
Dark Photons is quite difficult the LANL-P-25 PHENIX team is carefully planning
for the experiment. In Ref\cite{p25}, using EmCal hardware, one finds a
description of the hardware by Hubert van Hecke. Kun Liu, a co-PI, will plan
for the P-25 PHENIX team at Fermilab.

\section{Experimental Detection of Dark Matter}

Since we do not know what dark matter is, we need a diverse pool of
instruments and approaches to detect it\cite{Liu:2017drf,Cooley:2014aya}.
Experiments LUX~\cite{Akerib:2013tjd},
PandaX-II\cite{Tan:2016zwf}, PICO\cite{Amole:2016pye,Amole:2017dex},
DAMA/LIBRA\cite{Bernabei:2010mq,Bernabei:2013xsa,Bernabei:2018yyw},
SuperCDMS\cite{Agnese:2017njq}, CRESST-III~\cite{Abdelhameed:2019hmk},
which use direct detection methods, are reviewed.

The {\bf LUX} (Large Underground Experiment)~\cite{Akerib:2013tjd}
is a dark matter experiment that aims to detect the  WIMPs
(Weakly Interacting Massive Particles), a dark matter candidate. The
experiment is carried out 1.5 km underground under SURF (Sanford
Underground Research Facility). The LUX detector is a
two-phase xenon time-projection chamber (TPC), containing 370 kg of
xenon. The liquid xenon contains an array of
photosensors which can sense a single photon of light.

The second phase of the {\bf PandaX} (Particle and astrophysical Xenon)
project~\cite{Tan:2016zwf} is located in the China Jinping Underground
Laboratory (CJPL) and started it's operation in
early 2016. PandaX-II comprises of a 580 kg dual-phase xenon TPC along
with a 60x60 cm cylindrical target containing 55 top and 55 bottom
Hamamatsu R11410-20 3-inch photomultiplier tubes (PMTs).  The data sets
(Run 9 and Run 10)  are now used to perform
systematic studies using the PandaX-II detector.

The {\bf PICO} project\cite{Amole:2016pye} is working towards the detection of
dark matter particles with the bubble chamber technique. The expriments are
installed at a depth of 2 km in the SNOLAB underground laboratory at Sudbury,
Ontario, Canada. Searching WIMPs using superheated bubble chambers which are
operated in thermodynamic conditions at which they are almost insensitive to
gamma or beta radiation is the main detection principle of PICO experiment.
A bigger version of the experiment with up to 500 kg of active mass is in
progress.

The DAMA experiment at the Gran Sasso National Laboratories of the I.N.F.N. is
an observatory for rare processes that use radiopure scintillator set-ups.
The main activity field of {\bf DAMA/LIBRA}\cite{Bernabei:2010mq,
  Bernabei:2013xsa,Bernabei:2018yyw} is the investigation on Dark Matter
particles in the galactic halo.  The second generation DAMA/LIBRA set-up
(Large sodium Iodide Bulk for RAre processes is continuing
the investigations of DAMA/NaI.

The {\bf SuperCDMS}(SCDMS)\cite{Agnese:2017njq} is the
next generation of the CDMS(Cryogenic Dark Matter Search) II experiment,
which was located underground in the deep Soudan mine of Minnesota, USA.
SCDMS wants to
re-locate at SNOLAB (Vale Inco Mine, Sudbury, Canada) which is a much deeper
facility. CDMS-II comprises of a cryostat surrounded by a passive shield and
an outer muon veto situated beneath an overburden of 2090 meters water
equivalent. The use of underground facilities provide the required shielding
from the cosmic events and hence reduce the interference of known background
particles.

The {\bf CRESST}(Cryogenic Rare Event Search with Superconducting Thermometers)
experiment\cite{Abdelhameed:2019hmk}, located at the Gran Sasso underground
laboratory in Italy, has operated 13 detector modules in CRESST-III Phase 1
(2016-2018). CRESST searches for dark matter particles
via their elastic scattering off the nuclei in a target material.
The CRESST target comprises of scintillating $CaWO_{4}$ crystals.This
experiment searches for dark matter masses below 1.7 $GeV/c^2$, extending
the CRESST-II result in 2015. For CRESST-III, whose Phase 1 started in July
2016, detectors have been optomized to further probe the low-mass region.

The {\bf DINO} Dark Matter experiment\cite{dino16} will search for WIMPS
using scintillation detectors for the detection of recoiling nuclei.

The presence of low mass WIMPs in the various ongoing
experiments~\cite{lowmass1} have generated motivations to
cross-check future proposals\cite{Seth:2019sff}.
\newpage

\section{Conclusions}

In this Review we have reviewed many aspects of experiments and theory related
to Dark Matter.

Cosmic Microwave Background Radiation (CMBR) experiments have shown that about
23 \% of the Universe is Dark Matter and only about 4 \% of the universe is
normal matter.

The creation of Dark Matter mass during the EWPT and QCDPT, two Cosmological
Phase Transistions, via the dark matter field interacting with the dark
enargy (quintessence) field are reviewed. The Dark Matter mass $M_{DM}$ created
during the QCDPT could be more than 1 TeV,  more than an order of magnitude
larger than $M_{DM}$ created during the EWPT.

Theoretical and experimental studies of neutrino oscillations producing
sterile neutrinos are reviewed, with massive sterile neutrinos possibly
being the Dark Matter detected by CMBR experiments. Also a Dark Photon search,
with photons vector bosons rather than sterile neutrinos, which are fermions,
is briefly reviewed.

Finally, six experiments for the experimental detection of Dark Matter are
reviewed, with a detailed explanation of how these experiments can detect
Dark Matter.

\vspace{5mm}
\Large
{\bf Acknowledgements
\normalsize
\vspace{5mm}

Author D.D. acknowledges the facilities of Saha Institute of Nuclear Physics,
Kolkata, India. Author L.S.K. acknowledges support in part as a visitor at Los
Alamos National Laboratory, Group P25.
\vspace{5mm}