\subsection{Numerical methods}
\label{sec:kin_numerics}


Before we begin, the goal of the Section is clearly stated: We wish to evaluate the double integral
\begin{equation}
    W_{\ell, r} = 2 \int_0^\infty \int_0^\infty \exp{\lrp-\sg^2\rrp} \sg^{2r + 3} \lsp 1 - \cos^\ell\lrp\chi\lrp \sg, b\rrp\rrp\rsp \pfrac{b}{\sigma} \d \pfrac{b}{\sigma} \d \sg
    \label{eq:kin_dimless_collision_integral_method}
\end{equation}
numerically. To do so, we must evaluate the deflection angle $\chi(\sg, b; T)$, which is computed as $\chi(\sg, b; T) = \pi - \theta(R; \sg, b, T)$, where
\begin{equation}
    \begin{split}
        \theta(R; \sg, b, T) &= \int_R^{\infty} \lrp \frac{r^4}{b^2}\lsp 1 - \frac{u_{12}(r)}{k_B T \sg^2}\rsp - r^2 \rrp^{-\frac{1}{2}} \d r.
    \end{split}
    \label{eq:kin_dyn_theta_integral_method}
\end{equation}
The lower limit of integration, $R$, is the distance of closest approach, obtained by solving the Equation
\begin{equation}
         \frac{R^4}{b^2}\lsp 1 - \frac{u_{12}(R)}{k_B T \sg^2}\rsp - R^2 = 0.
    \label{eq:kin_dyn_R_root_method}
\end{equation}

The complete method therefore consists of several discrete parts: A root solver to determine $R$; a numerical integration scheme to determine $\theta(R; \sg, b, T)$, and thereby $\chi(\sg, b; T)$; and finally a 2d numerical integration scheme to evaluate $W_{\ell, r}$. Note that the temperature, $T$, is treated as a parametric dependency, meaning that $W_{\ell, r}$ must be computed at a given temperature.

These functions and equations are to a certain degree badly behaved, and applying commonly known integration schemes and root solvers was not directly successful. A part of the problem lies in the fact that the behaviour of $\frac{\d \theta}{\d r}(r; \sg, b, T)$ and Equation \eqref{eq:kin_dyn_R_root_method}, which must be solved to obtain a value for $R$, depend upon the molecular parameters, as well as having a parametric dependency on $\sg$, $b$ and $T$. Therefore robust algorithms for solving Equation \eqref{eq:kin_dyn_R_root_method} and evaluating the integral of Equation \eqref{eq:kin_dyn_theta_integral_method}, which is required to obtain a value for $\chi$, that are reliable for a large set of parameters are required. Further, the integrand of Equation \eqref{eq:kin_dimless_collision_integral_method}, the dimensionless collision integral, is costly to evaluate and badly behaved. This section describes the integration and root solving algorithms implemented to handle these issues.

\subsubsection{Solving for \mathsec{R}}
When implementing Newtons method to solve Equation \eqref{eq:kin_dyn_R_root_method} it is desirable to express the equation in such a way that the left hand side is either strictly increasing or decreasing in a wide region around the root. This makes the method less dependent upon a good initial guess, which is hard to obtain as the root is a function of the parameters $\sg$ and $b$, as well as the temperature. Analysing the Equation, one finds that this can be achieved by dividing by $\frac{R^4}{b^2}$ to obtain
\begin{equation}
    \psi_R(R; g, b, T) \equiv 1 - \frac{u_{12}(R)}{k_B T \sg^2} - \pfrac{b}{R}^2 = 0
    \label{eq:kin_num_R_root}
\end{equation}
Where $\psi_R$ is only a placeholder to refer to the function whose root determines $R$. Differentiation of $\psi_R$ with respect to $R$ yields
\begin{equation}
    \psi_R'(R; g, b, T) = - \frac{u'_{12}(R)}{k_B T \sg^2} + \frac{2b^2}{R^3}.
\end{equation}

For $b >> \sigma_{12}$, there will be very little interaction between the particles, except for at very small $T \sg$ ($T \sg < 150$). Therefore $b$ is used as an initial guess. However, note that due to the attractive part of the interaction potential, $\psi_R$ may have a local minima at $r < b$, with the root residing at a yet lower value of $r$. This would result in the Newton iteration continuing to infinity. Therefore, the implemented algorithm checks the second derivative of $\psi_R$ in each step, if the second derivative is negative, the iteration is restarted after reducing the initial guess by a factor of $0.1$. This was tested and found to be a robust and efficient method of finding the root of Equation \eqref{eq:kin_num_R_root}.

\subsubsection{Evaluating \mathsec{$\theta$}}
It is immediately clear that the integrand of Equation \eqref{eq:kin_dyn_theta_integral_method} is divergent at $r = R$, as this is implied by the solution of equation \eqref{eq:kin_dyn_R_root_method}. Further, as shown in Figure \ref{fig:theta_integration}, the integrand decreases rapidly at first, but does not converge very quickly. The algorithm used to evaluate the integral implements a trapezoidal integration rule, denoting $\psi_{\theta}(r; \sg, b, T) \equiv \frac{\d \theta}{\d r}$,
\begin{equation}
    I[r_a, r_b] = \sum_i \frac{\psi_{\theta}(r_i) + \psi_{\theta}(r_{i+1})}{2} \Delta r_i \approx \int_{r_a}^{r_b} \psi_{\theta} \d r.
\end{equation}
Where $I[r_a, r_b]$ denotes the numerical approximation to the integral on the interval from $r_a$ to $r_b$. This integration rule was chosen as it both gives solid control of the integration error and is fairly simple to implement. Due to the large second derivative of $\psi_{\theta}$ with respect to $r$ when $\vert r - R \vert \approx 0$, high resolution is required at small $r$. However, how high the resolution must be varies greatly with the parameters $\sg$, $b$ and $T$. Therefore it was made use of the fact that for the trapezoidal integration scheme
\begin{equation}
    \epsilon_{abs} = \left| I[r_a, r_b] - \int_{r_a}^{r_b} \psi_{\theta} \d r \right| \leq  \frac{\stackmax{(r_a, r_b)}\lcp \vert \psi_{\theta}''(r) \vert\rcp (r_b - r_a)^3}{12 N^2} 
\end{equation}
where $\epsilon_{abs}$ denotes the absolute error in the integral, $\stackmax{(r_a, r_b)}\lcp \vert \psi_{\theta}''(r) \vert\rcp$ denotes the maximum absolute value of the second derivative of the integrand on the interval $(r_a, r_b)$, and $N$ denotes the number of intervals. Utilising this, a tolerance, $\epsilon_{tol}^{step}$, may be selected for each interval, and the length of that interval can be determined as
\begin{equation}
    \Delta r_i = \lrp \frac{12 \epsilon_{tol}^{step}}{\stackmax{(r_i, r_i + \Delta r_i)}\lcp \vert \psi_{\theta}''(r)\vert\rcp}\rrp^{\frac{1}{3}}
\end{equation}
ensuring that the error over a single integration step never exceeds $\epsilon_{tol}^{step}$. The second derivative of the integrand of Equation \eqref{eq:kin_dyn_theta_integral_method}, i.e. $\frac{\d^3 \theta}{\d r^3}$ is obtained by straight forward differentiation as
\begin{equation}
    \begin{split}
        \psi_{\theta}''(r) &= \frac{3}{4}\varphi^{-\frac{5}{2}} \varphi'(r)- \frac{1}{2}\varphi^{-\frac{3}{2}} \varphi''(r)\\
        \varphi &= \frac{r^4}{b^2}\left(1 - \frac{u_{12}(r)}{k_B T \sg}\right) - r^2\\
        \varphi'(r) &= \frac{4r^3}{b^2} \left(1 - \frac{u_{12}(r)}{k_B T \sg}\right) + \frac{r^4}{b^2}\frac{u'_{12}(r)}{k_B T \sg} - 2r\\
        \varphi''(r) &= \frac{12 r^2}{b^2}\left(1 - \frac{u_{12}(r)}{k_B T \sg}\right) + \frac{8r^3}{b^2}\frac{u_{12}'(r)}{k_B T \sg} + \frac{r^4}{b^2}\frac{u_{12}''(r)}{k_B T \sg} -2\\
    \end{split}
\end{equation}
Noting that the second derivative is strictly decreasing in the most critical region (at small $r$), the step length is computed as
\begin{equation}
    \Delta r_i = \lrp \frac{12 \epsilon_{tol}^{step}}{\vert \psi_{\theta}''(r_a) \vert}\rrp^{\frac{1}{3}}.
\end{equation}

To assess whether the integral has converged at a given point, the integration algorithm moves forward iteratively. In the first iteration, 100 integration steps are conducted, as this was found to be slightly less that what was most often required to achieve convergence. In the following iterations, the value of the integral is stored and $N_{0.1}$ integration steps are conducted, where $N_{0.1} = \frac{N_t}{9}$ and $N_t$ is the total number of integration steps that have been conducted in previous iterations. This means that if $N_i$ denotes the total number of integration steps conducted after iteration $i$, $N_1 = 0.9N_2 = (0.9)^2 N_3 = ... = (0.9)^{i - 1}N_i$. The algorithm terminates when the change in the integral after an iteration is less than a set tolerance, $\epsilon_{conv}$. The integration is illustrated in Figure \ref{fig:theta_integration}.

\begin{figure}
    \centering
    \includegraphics[width=.7\textwidth]{method/figures/theta_integration.eps}
    \caption{Illustration of the need for high refinement of the integral of Equation \eqref{eq:kin_dyn_theta_integral_method}, and the manner in which convergence is determined.}
    \label{fig:theta_integration}
\end{figure}

After termination, the relative error of the total integral is checked, and if exceeding a tolerance $\epsilon_{tol}^{tot}$, the step error tolerance $\epsilon_{tol}^{step}$ is reduced, and the integral is recomputed.

\subsubsection{Evaluating \mathsec{$W_{\ell,r}$}}
\label{subsubsec:eval_W}
The integrand of Equation \eqref{eq:kin_dimless_collision_integral_method} is expensive to evaluate, and as shown in Figure \ref{fig:dimless_integral_hard} high resolution is required in certain regions to achieve satisfactory precision due to large first and second derivatives. It is unfeasible to conduct the integration with this high resolution over a region large enough to ensure that the integral has converged, without a higher capacity for parallelization than what was possible with the tools available during this work. Additionally, the shape of the integrand varies with the parameters $\ell$, $r$ and $T$ as well as the molecular parameters. The evaluation of a single diffusion coefficient, at a given temperature, to the first Enskog approximation, requires that the integral of Equation \eqref{eq:kin_dimless_collision_integral_method} is evaluated nine times for different combinations of $\ell, r$, and the number increases rapidly with increasing order of approximation. This makes it unfeasible to manually consider each integral. The implementation developed in this work is written to utilise two processor cores, with the option of using up to eight, to simultaneously evaluate several integrals. Additionally, a solution for storing the value of evaluated integrals for later use was implemented.

\begin{figure}
    \centering
    \includegraphics{method/figures/integrand.eps}
    \caption{Illustrative example of the form of the integrand of Equation \eqref{eq:kin_dimless_collision_integral_method}. Interpolating the integrand in such a manner that one captures the ''valley'' between the peaks is non-trivial, as the form of said valley is dependent on the potential parameters and the temperature, as well as the parameters $r$ and $\ell$. The function has been triangulated by the adaptive meshing scheme described in Section \ref{subsubsec:eval_W}, with exaggerated refinement to clearly illustrate the form of the function. Note that the surface consists of triangles, areas that look like parallelograms consist of two bordering triangles with the same colour.}
    \label{fig:dimless_integral_hard}
\end{figure}

To achieve reliable and efficient evaluation of the dimensionless collision integral an adaptive meshing algorithm was developed and implemented. This was used in combination with a planar linear interpolation scheme. 

\subsubsection*{Planar linear interpolation}
Given a set of three points, $P = \lcp p_1, p_2, p_3 \rcp \subset \mathbb{R}^3$, denote the corresponding set of points in the $(x,y)$-plane as $\Tilde{P} = \lcp\Tilde{p}_1, \Tilde{p}_2, \Tilde{p}_3\rcp \subset \mathbb{R}^2$. A plane can be constructed to interpolate the integrand in $P$, express this interpolating plane as 

\begin{equation}
    z(x, y) = Ax + By + C
\end{equation}

Further, ensure $x_1 \leq x_2 \leq x_3$, by reordering the points if necessary. Now, determine the functional form of the lines connecting each pair of points in $\Tilde{P}$, denoted $l_{ij}$ for the line passing through $(\Tilde{p}_i, \Tilde{p}_j)$. These have the form
\begin{equation}
    l_{ij} = A_{ij}^{(l)}x + B_{ij}^{(l)}, \hspace{.5cm} A_{ij}^{(l)} = \frac{y_j - y_i}{x_j - x_i}, \hspace{.5cm} B_{ij}^{(l)} = y_j - A_{ij}^{(l)} x_j
\end{equation}
The integral of the interpolating plane may now be written as
\begin{align}
    I[p_1, p_2, p_3] &= \int_{x_1}^{x_2}\int_{l_{12}}^{l_{13}} z(x,y) \d y \d x + \int_{x_2}^{x_3}\int_{l_{12}}^{l_{23}} z(x,y) \d y \d x, & y_2 &< y_3 \lor (y_2 = y_3 \land y_2 > y1)\\
    I[p_1, p_2, p_3] &= \int_{x_1}^{x_2}\int_{l_{13}}^{l_{12}} z(x,y) \d y \d x + \int_{x_2}^{x_3}\int_{l_{13}}^{l_{23}} z(x,y) \d y \d x, & y_2 &> y_3 \lor (y_2 = y_3 \land y_2 < y1)
\end{align}
Only the first term of the first case is treated explicitly here, the remaining terms are obtained by appropriate index-swapping. 
\begin{equation}
    \begin{split}
        \int_{x_1}^{x_2} \int_{l_{12}}^{l_{13}} z(x,y) \d y \d x &= \int_{x_1}^{x_2} \lsp \lrp Ax + C \rrp y + \frac{B}{2}y^2 \rsp_{l_{12}}^{l_{13}}\d x\\
        &= \int_{x_1}^{x_2} \lrp Ax + C \rrp \lrp \lrp A_{13}^{(l)} - A_{12}^{(l)}\rrp x + B_{13}^{(l)} - B_{12}^{(l)} \rrp \\
        &\hspace{1.5cm}+ \frac{B}{2}\lsp \lrp A_{13}^{(l)} x + B_{13}^{(l)} \rrp^2 - \lrp A_{12}^{(l)} x + B_{12}^{(l)} \rrp ^2 \rsp \d x\\
        & \equiv \int_{x_1}^{x_2} \Tilde{A}_{12}^{(13)}x^2 + \Tilde{B}_{12}^{(13)} x + \Tilde{C}_{12}^{(13)} \d x
    \end{split}
\end{equation}
where
\begin{equation}
    \begin{split}
        \Tilde{A}_{12}^{(13)} &= A \lrp A_{13}^{(l)} - A_{12}^{(l)}\rrp + \frac{B}{2} \lsp \lrp A_{13}^{(l)} \rrp^2 - \lrp A_{12}^{(l)} \rrp^2 \rsp\\
        \Tilde{B}_{12}^{(13)} &= A \lrp B_{13}^{(l)} - B_{12}^{(l)} \rrp + B \lrp A_{13}^{(l)} B_{13}^{(l)} -  A_{12}^{(l)} B_{12}^{(l)} \rrp + C \lrp A_{13}^{(l)} - A_{12}^{(l)}\rrp\\
        \Tilde{C}_{12}^{(13)} &= \frac{B}{2} \lsp \lrp B_{13}^{(l)} \rrp^2 - \lrp B_{12}^{(l)} \rrp^2 \rsp + C \lrp B_{13}^{(l)} - B_{12}^{(l)} \rrp
    \end{split}
\end{equation}
The integral in question is then
\begin{equation}
    \int_{x_1}^{x_2} \int_{l_{12}}^{l_{13}} z(x,y) \d y \d x = \frac{\Tilde{A}_{12}^{(13)}}{3}\lrp x_2^3 - x_1^3 \rrp + \frac{\Tilde{B}_{12}^{(13)}}{2} \lrp x_2^2 - x_1^2 \rrp + \Tilde{C}_{12}^{(13)} \lrp x_2 - x_1 \rrp.
\end{equation}

\subsubsection*{Adaptive meshing}
The implemented integration scheme is exact for any function $\psi(x, y) $ that satisfies $ \spder{\psi}{x}{x} = \spder{\psi}{x}{y} = \spder{\psi}{y}{y} = 0$ over the entire integration domain. It is assumed that, similarly to the analogous one-dimensional trapezoid scheme, the error in the two dimensional scheme scales linearly with the second derivative, and is proportional to the cube of the step length. That is, the absolute integration error in each subinterval $(x, x + h_x) \times (y, y+h_y)$ is assumed to be of the form

\begin{equation}
    \abs{\epsilon_i^{abs.}} \sim \abs{\spder{f}{x}{x} h_x^3} + \abs{\spder{f}{y}{y} h_y^3} + \abs{\spder{f}{x}{y} (h_x + h_y)^3}
\end{equation}

The following meshing algorithm increases the mesh resolution in regions with high second derivatives to ensure precise integration. At the same time, large step sizes are permitted in regions with small second derivatives to prevent unnecessary function evaluations from making the integration unnecessarily computationally intensive. 

Let $\mathcal{S} = \lcp (x_0, x_0 + \Delta x) \times (y_0, y_0 + \Delta y) \rcp$ be the total domain of integration. Further, let 
\begin{equation}
    \begin{split}
        \mathcal{S}^* &= \lcp (x_i^{\mathcal{S}^*}, y_i^{\mathcal{S}^*}, h_x^{\mathcal{S}^*}, h_y^{\mathcal{S}^*})\ |\ (x_i^{\mathcal{S}^*}, y_i^{\mathcal{S}^*}) \subseteq \mathcal{S}\ \land\ h_x^{\mathcal{S}^*} = x_{i + 1}^{\mathcal{S}^*} - x_i^{\mathcal{S}^*},\ h_y^{\mathcal{S}^*} = y_{j+1}^{\mathcal{S}^*} - y_j^{\mathcal{S}^*}\rcp\\
        \mathcal{S}_{ij}^*\lrp\mathcal{S}^*\rrp &= \lcp (x, y)\ \in \ \lrp x_i^{\mathcal{S}^*}, x_i^{\mathcal{S}^*} + h_x^{\mathcal{S}^*}\rrp \times \lrp y_j^{\mathcal{S}^*}, y_j^{\mathcal{S}^*} + h_y^{\mathcal{S}^*}\rrp\rcp\\
        \mathcal{S}_{ij}^{\Delta}(h_x', h_y', \mathcal{S}^*) &= \lcp (x_{i'}, y_{j'}) \subset \mathcal{S}_{ij}^*\lrp\mathcal{S}^*\rrp\ |\ x_{i' + 1} - x_{i'} = h_x'\ \land\ y_{j' + 1} - y_{j'} = h_y' \rcp
    \end{split}
    \label{eq:integration_domains}
\end{equation}
Such that $\mathcal{S}^*$ is some discrete subset of $\mathcal{S}$, $\mathcal{S}_{ij}^*$ is a continuous subset of $\mathcal{S}$, bounded by $\mathcal{S}^*$, and $\mathcal{S}_{ij}^{\Delta}$ is a discrete subset of $\mathcal{S}_{ij}^*$. The mappings
\begin{align}
    \mathcal{S}_{ij}^{x}\lrp\mathcal{S}^*\rrp\ :\ \mathcal{S}^* \mapsto &\left\{(x_k, y_{\ell}, h_x', h_y')\ |\ (x_{k}, y_{\ell}) \in \ \mathcal{S}_{ij}^{\Delta}(h_x', h_y', \mathcal{S}^*) \ :\ h_x' = \frac{h_x^{\mathcal{S}^*}}{2},\ h_y' = h_y^{\mathcal{S}^*} \right\}\\
    \mathcal{S}_{ij}^{y}\lrp\mathcal{S}^*\rrp\ :\ \mathcal{S}^* \mapsto &\left\{(x_k, y_{\ell}, h_x', h_y')\ |\ (x_{k}, y_{\ell}) \in \ \mathcal{S}_{ij}^{\Delta}(h_x', h_y', \mathcal{S}^*) \ :\ h_x' = h_x^{\mathcal{S}^*},\ h_y' = \frac{h_y^{\mathcal{S}^*}}{2} \right\}\\
    \mathcal{S}_{ij}^{xy}\lrp\mathcal{S}^*\rrp\ :\ \mathcal{S}^* \mapsto &\left\{(x_k, y_{\ell}, h_x', h_y')\ |\ (x_{k}, y_{\ell}) \in \ \mathcal{S}_{ij}^{\Delta}(h_x', h_y', \mathcal{S}^*) \ :\ h_x' = \frac{h_x^{\mathcal{S}^*}}{2},\ h_y' = \frac{h_y^{\mathcal{S}^*}}{2} \right\}
    \label{eq:integration_refine_maps}
\end{align}
may then be used to recursively generate discrete subsets of $\mathcal{S}$ with increased resolution in each iteration, the domains and mappings are illustrated in Figure \ref{fig:integration_domains}.
\begin{figure}[htb]
    \centering
    \includegraphics[width=\textwidth]{kinetic_gas/figures/integration_domains.eps}
    \caption{Domains and mappings of Equations \eqref{eq:integration_domains}-\eqref{eq:integration_refine_maps}, used to recursively refine the mesh used to interpolate the integrand of Equation \eqref{eq:kin_dimless_collision_integral_method}.}
    \label{fig:integration_domains}
\end{figure}

To integrate a function $\psi(x, y)$ over $\mathcal{S}$, the integration algorithm starts with a large step size in each direction, denoted $h_x = N_x \delta x$ and $h_y = N_y \delta y$, where $(N_x, N_y) = \lcp (2^{n_x}, 2^{n_y}) : (n_x, n_y) \in \mathbb{N}^2 \rcp$ and $\delta x, \delta y$ denote the minimum step size (i.e. maximum refinement) in each direction. The numerical second derivatives of $\psi$ are computed as 
\begin{equation}
    \begin{split}
        \spder{\psi}{x}{x} &\approx \lrp \Delta_{x}^2 \psi \rrp_{i,j} \equiv \frac{\psi_{i,j} - 2 \psi_{i+1, j} + \psi_{i + 2, j}}{(h_x)^2}\\
        \spder{\psi}{y}{y} &\approx \lrp \Delta_{y}^2 \psi \rrp_{i,j} \equiv \frac{\psi_{i,j} - 2 \psi_{i, j+1} + \psi_{i, j + 2}}{(h_y)^2}\\
        \spder{\psi}{x}{y} &\approx \lrp \Delta_{xy}^2 \psi\rrp_{i,j} \equiv \frac{\psi_{i, j} - 2 \psi_{i+1,j+1} + \psi_{i + 2, j + 2} - h_x^2 \Delta_x^2 \psi - h_y^2\Delta_y^2 \psi}{2 h_x h_y}
    \end{split}
\end{equation}

where $\psi_{k,l} = \psi(x_0 + kh_x, y_0 + lh_y)$ and the equivalencies define the $\Delta_{x}^2\psi$, $\Delta_y^2 \psi$ and $\Delta_{xy}^2\psi$ - finite difference operators. Once the numerical second derivatives are computed, the conditions of Equations \eqref{eq:refine_condition_1}-\eqref{eq:refine_condition_5} are checked successively.

\begin{align}
    &(\mathrm{I}) & \abs{\lrp \Delta_{x}^2 z \rrp_{i,j}} > \epsilon \hspace{.25cm}\land \hspace{.25cm}  \abs{\lrp \Delta_{y}^2 z \rrp_{i,j}} < \epsilon  \hspace{.25cm} \land  \hspace{.25cm} \abs{\lrp \Delta_{xy}^2 z\rrp_{i,j}}  < \epsilon \hspace{.25cm}\land \hspace{.25cm} N_x > 1, \label{eq:refine_condition_1}\\
    &(\mathrm{II})& \abs{\lrp \Delta_{x}^2 z \rrp_{i,j}} < \epsilon \hspace{.25cm}\land \hspace{.25cm}  \abs{\lrp \Delta_{y}^2 z \rrp_{i,j}} > \epsilon  \hspace{.25cm} \land  \hspace{.25cm} \abs{\lrp \Delta_{xy}^2 z\rrp_{i,j}}  < \epsilon \hspace{.25cm}\land \hspace{.25cm} N_y > 1,\\
    &(\mathrm{III})& \abs{\lrp \Delta_{x}^2 z \rrp_{i,j}} +  \abs{\lrp \Delta_{y}^2 z \rrp_{i,j}} + \abs{\lrp \Delta_{xy}^2 z\rrp_{i,j}}  > \epsilon \hspace{.25cm}\land \hspace{.25cm} N_x > 1 \hspace{.25cm}\land \hspace{.25cm} N_y > 1,\\
    &(\mathrm{IV})& \abs{\lrp \Delta_{x}^2 z \rrp_{i,j}} +  \abs{\lrp \Delta_{y}^2 z \rrp_{i,j}} + \abs{\lrp \Delta_{xy}^2 z\rrp_{i,j}}  > \epsilon \hspace{.25cm}\land \hspace{.25cm} N_x > 1,\\
    &(\mathrm{V})& \abs{\lrp \Delta_{x}^2 z \rrp_{i,j}} +  \abs{\lrp \Delta_{y}^2 z \rrp_{i,j}} + \abs{\lrp \Delta_{xy}^2 z\rrp_{i,j}}  > \epsilon \hspace{.25cm}\land \hspace{.25cm} N_y > 1, \label{eq:refine_condition_5}
\end{align}

Here, $\epsilon$ is some predetermined tolerance value. Upon reaching a condition that is true, the algorithm selects a map, following the rules
\begin{align*}
    (\mathrm{I}) \hspace{.15cm} \lor (\mathrm{IV}) \mapsto  \mathcal{S}_{ij}^{x},& &(\mathrm{II}) \hspace{.15cm} \lor (\mathrm{V}) \mapsto \mathcal{S}_{ij}^{y},& &(\mathrm{III}) \mapsto \mathcal{S}_{ij}^{xy},\\
\end{align*}
to generate a subdomain and calls itself recursively on that subdomain.

If no condition is true, the algorithm steps across the current domain, following the pattern $p_0 = (x_0, y_0)$, $p_1 = (x_0, y_1)$, $p_2 = (x_1, y_0)$, ..., $p_{2k} = (x_{k}, y_0)$, $p_{2k + 1} = (x_{k}, y_1)$. Upon reaching the edge of the current domain, at $x_k = x_N$, two steps are conducted in the $y$-direction, and the step in the $x$-direction is reversed, such that the stepping follows the pattern $p_{N + 2k} = (x_{N - k}, y_1)$, $p_{N + 2k + 1} = (x_{N - k}, y_2)$. The pattern is illustrated in Figure \ref{fig:stepping_pattern}.
\begin{figure}
    \centering
    \includegraphics{kinetic_gas/figures/stepping_algorithm.eps}
    \caption{The stepping algorithm used to generate the recursively refined integration mesh. Arrows indicate integration steps, numbers indicate the order in which the interpolating triangles are generated, going in the order $1a-12a-1b-12b-1c-20c$. Refined regions are indicated, with yellow regions refined by the map $\mathcal{S}_{ij}^{x}$, orange regions refined by $\mathcal{S}_{ij}^{y}$ and purple regions refined by $\mathcal{S}_{ij}^{xy}$. The region $3c-6c$ is doubly refined, first by $\mathcal{S}_{ij}^{x}$, then by $\mathcal{S}_{ij}^{y}$.}
    \label{fig:stepping_pattern}
\end{figure}

To verify the reliability of the implemented integration algorithm, several functions with analytic integrals were meshed and integrated numerically. As shown in Figure \ref{fig:demo_adaptive_mesh2}, the meshing algorithm successfully increased resolution in areas with large second derivatives. The resulting numeric integrals had errors of typically less than $3\ \%$, even with large initial grid size. To minimise the number of function evaluations, the integrands' value in each gridpoint is stored, such that it will not be recomputed for each triangle that uses the point.
\begin{figure}[htb]
    \centering
    \subfloat[$\epsilon_r = 0.38$ \%]{\includegraphics[width=.48\textwidth]{kinetic_gas/figures/adaptive_mesh_demo2.png}}
    \subfloat[$\epsilon_r = 0.48$ \%]{\includegraphics[width=.48\textwidth]{kinetic_gas/figures/adaptive_mesh_demo3.png}}
    \caption{Visualisation of the triangulation of two test functions. $\epsilon_r = \frac{\abs{I_{\text{exact}} - I_{\text{numeric}}}}{I_{\text{exact}}}$ is the relative error in the numeric integrals. In (a), notice the individual refinement in the x- and y- directions.}
    \label{fig:demo_adaptive_mesh2}
\end{figure}

\subsubsection*{Refining the algorithm parameters}
From geometric considerations one will find that for a hard sphere

\begin{equation}
    \chi^{HS} = \cos^{-1}\lsp 1 - 2 \lrp 1 - \pfrac{b}{\sigma_{ij}^{HS}}^2\rrp \rsp \times \begin{cases}
    1, & b < \sigma_{ij}^{HS}\\
    0, & b > \sigma_{ij}^{HS}
    \end{cases}
\end{equation}

After testing the implemented algorithm for the functions shown in Figure \ref{fig:demo_adaptive_mesh2}, appropriate parameters required to precisely and efficiently evaluate the dimensionless collision integral were determined. This was done by using the pseudo-hard-sphere potential,
\begin{equation}
    u_{p-HS}(r_{ij}) = \begin{cases}
    \pfrac{\sigma_{ij}}{r_{ij}}^{\varphi} + \frac{\varphi(\varphi + 1)}{2}\pfrac{r_{ij}}{\sigma_{ij}}^2 - \varphi \lrp\varphi - \frac{\varphi}{\sigma_{ij}} + 1\rrp r_{ij} & \\
    \hspace{1.5cm}  + \varphi\lrp (\varphi + 1) \sigma_{ij} - \frac{\varphi}{2} - \frac{3}{2}\rrp - 1 & r_{ij}<\sigma_{ij}\\
    0 & r_{ij} > \sigma_{ij}
    \end{cases}
    \label{eq:pseudo_hs_potential}
\end{equation}
to compute $\chi$, $R$ and $W_{\ell,r}$ numerically, and comparing the result to the known analytical solutions for hard spheres. This potential was chosen due to the root solver and integration scheme requiring continuous first and second derivatives to behave as intended, and the given potential satisfies
\begin{equation}
    u_{p-HS}(\sigma) = u'_{p-HS}(\sigma) = u''_{p-HS}(\sigma) = 0.
\end{equation}
An exponent of $\varphi = 20$ was chosen arbitrarily as a large exponent to force the potential to rise steeply. Using the integration parameters shown in Table \ref{tab:integration_params}, the integration error when compared to the analytic solutions for a hard sphere were consistently below $1$ \% for $\ell, r < 9$. 

\begin{table}[htb]
    \centering
    \caption{Parameters used with the integration scheme presented in Section \ref{sec:kin_numerics} throughout this work. $\epsilon_0$ refers to the tolerance for refinement, used in Equations \eqref{eq:refine_condition_1}-\eqref{eq:refine_condition_5}, $\delta$ is the minimum step size along each axis and $N$ is the initial step size in multiples of $\delta$. $\sg_N$ and $\pfrac{b}{\sigma}_N$ are the cutoff values of the integration domain.}
    \begin{tabular}{c c | c c}
    \toprule
        Parameter & Value & Parameter & value \\
    \midrule
        $N_{\sg}$ & 4 & $N_{\pfrac{b}{\sigma}}$ & 16 \\
        $\delta_{\sg}$ & 0.5 & $\delta_{\pfrac{b}{\sigma}}$ & 0.05 \\
        $\sg_N$ & 7.5 & $\pfrac{b}{\sigma}_N$ & 5\\
        $\epsilon_0$ & $10^{-5}$ & & \\
    \bottomrule
    \end{tabular}
    \label{tab:integration_params}
\end{table}

The entirety of the implementation developed in this work will be made available in a public repository on \href{https://github.com/vegardjervell/Kineticgas/}{GitHub}, and can be supplied upon request.\cite{jervell_git}