\section{\label{sec:intro}Introduction}
The study of rarefied gases covers numerous applications, including flows caused by evaporation and condensation, upper-atmospheric dynamics, modeling of airborne particles, the reflective and reactive properties of gases interacting with solid and liquid surfaces and so on. Rarefied gas flows are characterized by a dimensionless parameter, the Knudsen number ($\mathrm{Kn}$), which is the ratio of the mean free path $\lambda$ of the gas and a characteristic length scale $L$ in the problem.
For very small values of the Knudsen number ($\mathrm{Kn}\lesssim 0.01$), the classical continuum theories, namely the Euler and Navier--Stokes--Fourier (NSF) equations, are quite effective in capturing rarefaction effects but they fall short of doing so when the Knudsen number is not very small.
Although the NSF equations fail to capture several non-equilibrium phenomena (like non-homogeneity in pressure profile and unusual temperature dip in the Poiseuille flow \citep{OSA1989,TTS2009,RRPM2016}, heat flux direction opposite to the temperature gradient or the cross effects where heat flows from a low-temperature region to a high-temperature region \citep{ST2008,MRS2015,RTS2013}), yet, by exploiting the coupling among the thermodynamic forces and fluxes to form a closed system, the range of applicability of the new equations is enhanced in comparison to the NSF equations.
A model that exploits the coupling among the thermodynamic forces and fluxes to yield an improved set of constitutive relations for the stress and heat flux appearing in the conservation laws has recently been propounded by Rana et al.~[\onlinecite{RGS2018}].
The constitutive relations for the stress and heat flux obtained in this model are coupled through a coupling coefficient; hence they are referred to as the coupled constitutive relations (CCR), and the model wherefore is referred to as the CCR model.
In the linearized and steady state, the CCR model reduces to the linearized Grad 13-moment (G13) equations \cite{Grad1949} in the steady state as a special case, and on taking the coupling coefficient as zero, the CCR model reduces to the original NSF equations.
Owing to its simplicity and viable features, the CCR model has been applied successfully to some problems pertaining to rarefied gas flows \cite{RSCLS2021, MS2022}.
Although, there do exist other models, such as the regularized 13-moment (R13) equations \cite{ST2003, Struchtrup2005}, the regularized 26-moment (R26) equations \cite{GuEmerson2009}, etc., that can describe rarefied gas flows somewhat more accurately than the CCR model, especially for flows at moderate Knudsen numbers, we shall use the CCR model due to its simplicity in this work.

In this paper, we shall focus our attention to exploring rarefied gas flow problems in two dimensions (2D).
The reason for this is twofold:
firstly, for a symmetric uniform flow in three dimensions (3D), it is sufficient to study the problem in 2D, thanks to symmetry of the problem,
and secondly, there are some intriguing problems in 2D that do not arise in 3D, such as Stokes' paradox \citep{Lamb1932}, which states the non-existence of the steady-state solution to Stokes' equation in 2D.
Furthermore, we shall investigate the problems numerically using a truly meshless numerical technique introduced by Kupradze and Aleksidze \citep{KA1964} known as the method of fundamental solutions (MFS).

The MFS is a boundary-type meshfree approach in which an approximate solution of a (linear) boundary value problem is expressed as a linear combination of numerous singular functions, referred to as the fundamental solutions, and the
boundary conditions are satisfied at several locations on the boundary, referred to as the boundary nodes or collocation points, aiming to determine the unknown coefficients in the linear combination.
Apart from being time efficient due to reduced spatial dimension in boundary discretization, the quality of being free from integrals makes the MFS peerless among other meshfree methods (such as the boundary element method \citep{brebbia1994BEM}, finite point method \citep{finite_pt_method}, diffuse element method \citep{diffuse_element}, element-free Galerkin method \citep{element_free_galerkin}) that involve complex integrals.
The MFS has proven to be an efficient executable numerical scheme in various areas such as thermoelasticity, electromagnetics, electrostatics, wave scattering, inverse problems and fluid flow problems; see, e.g.~Refs.~[\onlinecite{LFS2021, BK2001, FKM2003, KLM2011, YJFMT2006, LC2016}].
Moreover, the MFS is also suitable for the analysis of problems involving shape optimization, moving boundary and unknown boundary \citep{FR2002, YJFMT2006, chen2008, SLL2012, Alves2009}, since the problems of modeling and satisfying boundary conditions are relatively simple for them.
All these advantages make the application of the MFS to the CCR model evidently favorable.

Several researchers have employed the MFS to solve the Helmholtz-, harmonic- and biharmonic-type boundary value problems in 2D as well as in 3D, see e.g.~Refs.~[\onlinecite{PKG1998,ML2005}].
The MFS works as a good numerical strategy if the fundamental solutions to the problem are predefined.
In the past few years, there has been a surge of interest in employing the MFS to various models for rarefied gas flows, for instance to the NSF, G13, R13 and CCR models~\citep{LC2016, CSRSL2017, RSCLS2021}, because the predefined fundamental solutions of the well-known equations, such as the Laplace, Helmholtz and biharmonic equations, can be exploited to determine the fundamental solutions for the NSF, G13, R13 and CCR models.
Nevertheless, all the works on the MFS for rarefied gas flows have investigated the problems in 3D only.
It is, however, important to note that the two-dimensional fundamental solutions for a model cannot be deduced directly from its three-dimensional counterpart due to the fact that the associated Green's functions are entirely different in 2D and 3D.
Therefore the main objectives of the paper are (i) to determine the fundamental solutions of the linearized CCR model in 2D
and (ii) to implement the determined fundamental solutions in a numerical framework.
To gauge the accuracy of the developed numerical framework, the obtained numerical results are also validated against those obtained with other models for a few problems existing in the literature.
We pick two internal-flow problems from Refs.~[\onlinecite{Onishi1977, Aoki_Sone_1989}] that have been investigated in these references with the linearized Bhatnagar--Gross--Krook (BGK) model \citep{BGK1954} [also referred to as the Boltzmann--Krook--Welander (BKW) kinetic model by some authors \cite{Onishi1977, Aoki_Sone_1989, Sone2007}].
In the first problem, the evaporation and condensation of a mildly rarefied vapor confined between two coaxial cylinders is studied while in the second problem,  a temperature-driven rarefied gas flow between two non-coaxial cylinders is investigated.

A well-known shortcoming of the MFS is that it yields the results that are highly sensitive toward the location of the singularity points (also referred to as the source points), and that a high accuracy is accompanied with an ill-conditioned collocation matrix \citep{Alves2009, CKL2016, WLQ2018, CH2020}.
Therefore, we also investigate the optimal location of singularities using an approach based on the effective condition number as discussed in Refs.~[\onlinecite{DML2009,
CNYC2023, WL2011}].

The remainder of the paper is structured as follows.
The linearized CCR model and the generalized boundary conditions associated with it are outlined in Sec.\,\ref{sec:balance_eqs}.
The fundamental solutions for the CCR model in 2D are determined in Sec.\,\ref{sec:solutions}.
The technique to apply the MFS by forming a system of equations for any arbitrary geometry is discussed in Sec.\,\ref{sec:system}.
The implementation of the MFS along with its validation (i) for the problem of a vapor flow between two coaxial cylinders is  demonstrated in Sec.\,\ref{sec:prob1}
and (ii) for the problem of a temperature-driven rarefied gas flow between two non-coaxial cylinders is discussed in Sec.\,\ref{sec:prob2}.
The location of singularities based on the effective condition number approach is examined in Sec.\,\ref{sec:location}.
The paper ends with conclusions and outlook in Sec.\,\ref{sec:conclusion}.
\section{\label{sec:balance_eqs}The linearized CCR model and boundary conditions}
The CCR model consists of the conservation laws---the balance equations for the mass, momentum and energy---closed with the constitutive relations for the stress and heat flux, which are coupled with each other through a coupling coefficient.
The full details of the CCR model can be found in Ref.~[\onlinecite{RGS2018}].
In this work, we require them in the linearized form.
To this end, we convert the mass, momentum and energy balance equations and the coupled constitutive relations into a linear-dimensionless form by assuming small perturbations in flow fields from their respective equilibrium values.
The velocity, stress and heat flux in the equilibrium state vanish whereas the density and temperature in the equilibrium state are constants $\tilde{\rho}_0$ and $\tilde{T}_0$, respectively.
The dimensionless perturbations in the density $\tilde{\rho}$ and temperature $\tilde{T}$ from their values in the equilibrium are given by
\begin{align}
\label{scaledrhoT}
\rho=\frac{\tilde{\rho}-\tilde{\rho}_0}{\tilde{\rho}_0}
\quad\text{and}\quad
T=\frac{\tilde{T}-\tilde{T}_0}{\tilde{T}_0},
\end{align}
respectively.
Similarly, the dimensionless perturbations in the velocity $\tilde{\bm{v}}$, stress tensor $\tilde{\bm{\sigma}}$ and heat flux $\tilde{\bm{q}}$ from their values in the equilibrium are given by
\begin{align}
\label{scalednoneqbquan}
\bm{v}=\frac{\tilde{\bm{v}}}{\sqrt{\tilde{\theta}_0}},
\quad
\bm{\sigma}=\frac{\tilde{\bm{\sigma}}}{\tilde{\rho}_0 \tilde{\theta}_0}
\quad\text{and}\quad
\bm{q}=\frac{\tilde{\bm{q}}}{\tilde{\rho}_0 (\tilde{\theta}_0)^{3/2}},
\end{align}
respectively, where $\tilde{\theta}_0 = \tilde{R}\tilde{T_0}$ with $\tilde{R}$ being the gas constant.
The linearized equation of state $p\approx\rho+T$ gives the dimensionless perturbation in the pressure from its equilibrium value $\tilde{p}_0=\tilde{\rho}_0\tilde{\theta}_0$.
For the sake of simplicity, the field variables with tilde are the quantities with dimensions while those without tilde are the dimensionless quantities throughout the paper.

Considering $\tilde{L}$ to be the characteristic length scale, the dimensionless position vector is $\bm{r}=\tilde{\bm{r}}/\tilde{L}$.
Inserting these dimensionless variables into the CCR model \cite{RGS2018} and dropping all nonlinear terms in the perturbed variables, one readily obtains the linear-dimensionless CCR model.
Here, we present them directly.
The linear-dimensionless mass, momentum and energy balance equations in the steady state read
\begin{align}
\label{mass_bal}
\bm{\nabla}\cdot\bm{v}=0,
\\
\label{mom_bal}
\bm{\nabla} p+\bm{\nabla}\cdot\bm{\sigma}=\bm{0},
\\
\label{energy_bal}
\bm{\nabla}\cdot\bm{q}=0,
\end{align}
and, to close the system \eqref{mass_bal}--\eqref{energy_bal}, the linearized coupled constitutive relations\citep{RGS2018}
\begin{align}
\label{ccr_1}
\bm{\sigma}&=-2\mathrm{Kn}\overline{\bm{\nabla} \bm{v}}-2\alpha_0 \mathrm{Kn}\overline{\bm{\nabla} \bm{q}},
\\
\label{ccr_2}
\bm{q}&=-\frac{c_p \mathrm{Kn}}{\mathrm{Pr}}(\bm{\nabla} T+\alpha_0\bm{\nabla}\cdot\bm{\sigma}),
\end{align}
where $\alpha_0$ is the coupling coefficient through which constitutive relations \eqref{ccr_1} and \eqref{ccr_2} are coupled;
$c_p = \tilde{c}_p / \tilde{R}$ with $\tilde{c}_p$ being specific heat capacity of the gas at a constant pressure; and
\begin{align}
\mathrm{Pr}=\frac{5\tilde{R}\tilde{\mu}_0}{2\tilde{\kappa}_0}
\quad\text{and}\quad
\mathrm{Kn}=\frac{\tilde{\mu}_0}{\tilde{\rho}_0\sqrt{\tilde{\theta}_0}\tilde{L}}
\end{align}
are the Prandtl number and Knudsen number, respectively with $\tilde{\mu}_0$ and $\tilde{\kappa}_0$ being the viscosity and thermal conductivity at the equilibrium state.
The quantities $\overline{\bm{\nabla} \bm{v}}$ and $\overline{\bm{\nabla} \bm{q}}$ in \eqref{ccr_1} are the symmetric-tracefree parts of the tensors $\bm{\nabla} \bm{v}$ and $\bm{\nabla} \bm{q}$, respectively.
For a vector $\bm{\psi}$, the symmetric-tracefree part of the tensor $\bm{\nabla} \bm{\psi}$ is defined as \cite{Gupta2020}
\begin{align}
\overline{\bm{\nabla} \bm{\psi}}&=\frac{1}{2}\Big[{\bm{\nabla} \bm{\psi}}+(\bm{\nabla} \bm{\psi})^\mathsf{T}\Big]-\frac{1}{d}(\bm{\nabla}\cdot\bm{\psi})\bm{I},
\end{align}
where $d$ is the dimension and $\bm{I}$ is the identity tensor in $d$-dimensions.
Since we shall be dealing with the problems in 2D, $d=2$ is fixed throughout the paper.
Furthermore, the dimensionless specific heat of a gas at a constant pressure is $c_p = (5+\mathsf{n})/2$, where $\mathsf{n}$ is a positive integer that accounts for the rotational degrees of freedom in a polyatomic gas.
For monatomic gases, there is no rotational degree of freedom; consequently, $\mathsf{n}=0$ and $c_p=5/2$ for monatomic gases.
We shall only deal with monatomic gases in this paper, and hence $c_p=5/2$ throughout this paper.
Equations \eqref{mass_bal}--\eqref{energy_bal} closed with \eqref{ccr_1} and \eqref{ccr_2} are referred to as the linearized CCR model.
For $\alpha_0=0$, the linearized CCR model reduces to the linearized NSF equations and for $\alpha_0=2/5$, the linearized CCR model reduces to the linearized G13 equations.
Since we shall be comparing the results obtained in the present work with those obtained with the BGK model, for which the Prandtl number is unity \cite{Struchtrup2005}, $\mathrm{Pr}=1$  throughout this paper.
Also, the parameter $\alpha_0$ is taken as $0.3197$, the value of $\alpha_0$ for hard sphere molecules \cite{RGS2018}, throughout this paper.

The thermodynamically-consistent boundary conditions complementing the linear CCR model have been derived in Ref.~[\onlinecite{RSCLS2021}].
For a problem in 3D, the boundary conditions complementing the linear CCR model are given in Eqs.~$(4.2a)$, $(4.2b)$, $(4.3a)$ and $(4.3b)$ of Ref.~[\onlinecite{RSCLS2021}].
Eqs.~$(4.2a)$ and $(4.2b)$ of Ref.~[\onlinecite{RSCLS2021}] are the boundary conditions on the normal components of the mass and heat fluxes, respectively, while Eqs.~$(4.3a)$ and $(4.3b)$ of Ref.~[\onlinecite{RSCLS2021}] are the boundary conditions on the shear stress---two conditions due to two tangential directions in 3D.
Since in 2D, there will be only one tangential direction, boundary condition $(4.3b)$ of Ref.~[\onlinecite{RSCLS2021}] is irrelevant in the present work and the superscript `$(1)$' can be dropped from the unit tangent vector $\bm{t}^{(1)}$ in $(4.3a)$ of Ref.~[\onlinecite{RSCLS2021}] for simplicity.
Thus, the linear-dimensionless boundary conditions for the linearized CCR model in 2D read \cite{RSCLS2021}
\begin{align}
\label{bc_1}
(\bm{v}-\bm{v}^I)\cdot \bm{n}=& -\eta_{11}(p-p_{\mathrm{sat}} + \bm{n}\cdot\bm{\sigma}\cdot\bm{n})
\nonumber\\
&+\eta_{12}(T-T^I+\alpha_0 \bm{n}\cdot\bm{\sigma}\cdot\bm{n}),
\end{align}
\begin{align}
\label{bc_2}
\bm{q}\cdot \bm{n} =& \medspace\eta_{12}(p-p_{\mathrm{sat}} + \bm{n}\cdot\bm{\sigma}\cdot\bm{n})
\nonumber\\
&-(\eta_{22}+2\tau_0)(T-T^I+\alpha_0 \bm{n}\cdot\bm{\sigma}\cdot\bm{n}),
\end{align}
\begin{align}
\label{bc_3}
\bm{t}\cdot\bm{\sigma}\cdot\bm{n} =& -\varsigma (\bm{v}-\bm{v}^I+\alpha_0 \bm{q})\cdot\bm{t},
\end{align}
where $\bm{n}$ and $\bm{t}$ are the unit normal and tangent vectors, respectively.
In boundary conditions \eqref{bc_1}--\eqref{bc_3}, $\eta_{ij}$'s,  for $i,j\in\{1,2\}$  are the Onsager reciprocity coefficients, which from Sone's asymptotic kinetic theory \citep{Sone2007} turn out to be
\begin{align}
\left.
\begin{aligned}
\label{coff}
\eta_{11}&=0.9134\sqrt{\frac{2}{\pi}}\frac{\vartheta}{2-\vartheta},
\\
\eta_{12}&=0.3915\sqrt{\frac{2}{\pi}}\frac{\vartheta}{2-\vartheta},
\\
\eta_{22}&=0.1678\sqrt{\frac{2}{\pi}}\frac{\vartheta}{2-\vartheta}
\end{aligned}
\right\}
\end{align}
under the assumption of the accommodation coefficient being unity (which also holds true for the diffuse reflection boundary condition).
The parameter $\vartheta$ in the above coefficients is the evaporation/condensation coefficient.
For canonical boundaries and phase-change boundaries, $\vartheta=0$ and $1$, respectively, are the largely accepted values of $\vartheta$ in the literature.
The temperature-jump and velocity-slip coefficients are given by \cite{RSCLS2021}
\begin{align}
\tau_0=0.8503\sqrt{\frac{2}{\pi}}\quad\text{and}\quad
\varsigma = 0.8798\sqrt{\frac{2}{\pi}},
\end{align}
respectively.
Furthermore,  $\bm{v}^I$, $T^I$ and $p_{\mathrm{sat}}$ in boundary conditions \eqref{bc_1}--\eqref{bc_3} represent the velocity, temperature and saturation pressure at the interface.
It is important to note that the coefficients $\alpha_0$ in boundary conditions \eqref{bc_1}--\eqref{bc_3} are actually the fitting parameters and could be different from the coupling coefficient $\alpha_0$.
Moreover, the coefficient $\alpha_0$ in each of boundary conditions \eqref{bc_1}--\eqref{bc_3} could also be different from each other.
The only reason that the coefficients $\alpha_0$ in boundary conditions \eqref{bc_1}--\eqref{bc_3} have been taken as the same as the coupling coefficient in the CCR model because the boundary conditions obtained in this way are thermodynamically consistent \cite{RGS2018}.

\section{\label{sec:solutions} Derivation of the fundamental solutions of the CCR model}
The fundamental solutions of the CCR model in 3D have already been derived in Ref.~[\onlinecite{RSCLS2021}].
However,  as mentioned in Sec.\,\ref{sec:intro},  the fundamental solutions of a model in 2D and 3D are independent of each other because the inherent Green's functions are independent of each other in 2D and 3D;
consequently,   the fundamental solution of a model in 2D cannot be determined from its 3D counterpart in general.
Therefore,  we derive the  fundamental solutions of the CCR model in 2D from scratch in this section.
To this end,  we add a Dirac delta forcing term of strength $\bm{f}$ on the right-hand side of the momentum balance equation to represent a point force, and a point heat source of strength $g$  on the right-hand side of the energy balance equation.
Furthermore,  to deal with phase-change effects at the liquid-vapor interface,  a point mass source of strength $h$ is also added  on the right-hand side of the mass balance equation.
For determining the fundamental solutions of a system of partial differential equations, it is customary to consider only one point source  at a time and then to superimpose the solutions obtained by taking each point source at a time in order to incorporate the effects of all point sources; see,  e.g. ~Refs.~[\onlinecite{LC2016, RSCLS2021}].
Nevertheless,  we take all three point sources $\bm{f}$,  $g$ and $h$ simultaneously and solve the resulting system of equations altogether.
We have verified---although not shown here for brevity---that this procedure also yields exactly the same solution as that obtained  by superimposing the solutions obtained by solving the systems separately with one point source at a time.
It is easy to keep track of the derivation of the fundamental solutions of the CCR equations in the indicial notation.
Therefore,  we shall derive the fundamental solutions of the CCR equations first in the indicial notation and then express them in the vectorial/tensorial notation.
For that,  let us first write down the linearized CCR model (with the point source terms) in the indicial notation.
The mass,  momentum and energy balance equations \eqref{mass_bal}--\eqref{energy_bal}---with the point source terms---in the indicial notation read
\begingroup
\allowdisplaybreaks
\begin{align}
\label{mass_bl}
\frac{\partial v_i}{\partial x_i}&=h\,\delta(\bm{r}),
\\
\label{mom_bl}
\frac{\partial p}{\partial x_i}+\frac{\partial \sigma_{ij}}{\partial x_j}&=f_i \,\delta(\bm{r}),
\\
\label{energy_bl}
\frac{\partial q_i}{\partial x_i}&=g\,\delta(\bm{r}),
\end{align}
\endgroup
where the Einstein summation is assumed over repeated indices in a term.
The CCR closure (Eqs.~\eqref{ccr_1} and \eqref{ccr_2}) in the indicial notation reads
\begin{align}
\label{ccr__1}
\sigma_{ij}=&-\mathrm{Kn}\left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}-\frac{\partial v_k}{\partial x_k}\delta_{ij}\right)\nonumber\\&
-\alpha_0\mathrm{Kn}\left(\frac{\partial q_i}{\partial x_j}+\frac{\partial q_j}{\partial x_i}-\frac{\partial q_k}{\partial x_k}\delta_{ij}\right),
\end{align}
\begin{align}
\label{ccr__2}
q_i=&-\frac{c_p \mathrm{Kn}}{\mathrm{Pr}}\left(\frac{\partial T}{\partial x_i}+\alpha_0 \frac{\partial \sigma_{ij}}{\partial x_j}\right).
\end{align}
Note that $\delta_{ij}$ in Eq.~\eqref{ccr__1} is the Kronecker delta and that $d=2$ has been substituted while writing Eq.~\eqref{ccr__1}.
We solve system \eqref{mass_bl}--\eqref{ccr__2} using the Fourier transform, which for a function $\hat{F}(\bm{\omega})$ is defined as
\begin{align}
\label{ft}
\mathcal{F}\big(F(\bm{r})\big)
=\hat{F}(\bm{\omega}) = \int_{\mathbb{R}^2} F(\bm{r}) \, \mathrm{e}^{\mathbbm{i} \, \bm{\omega} \cdot \bm{r}} \, \mathrm{d} \bm{r}
\end{align}
and the corresponding inverse Fourier transformation is defined as
\begin{align}
\label{ftinv}
\mathcal{F}^{-1}\big(\hat{F}(\bm{\omega})\big)
=F(\bm{r})
=\frac{1}{(2\pi)^2}\displaystyle{ \int_{\mathbb{R}^2}}
\hat{F}(\bm{\omega}) \, \mathrm{e}^{-\mathbbm{i}\,\bm{\omega} \cdot \bm{r}} \, \mathrm{d}\bm{\omega}.
\end{align}
Applying the Fourier transformation in Eqs.~\eqref{mass_bl}--\eqref{ccr__2} and using the fact that $\mathcal{F}[\delta(\bm{r})]=1$,  we obtain
\begin{align}
\label{mass_ft}
 \omega_i \hat{v}_i &=\mathbf{\mathbbm{i}} h,
 \\
 \label{mom_ft}
 \omega_i \hat{p}+\omega_j\hat{\sigma}_{ij} &= \mathbbm{i} f_i,
\\
\label{energy_ft}
 \omega_i \hat{q}_i &= \mathbbm{i} g,
 \end{align}
\begin{align}
 \label{ccr1_ft}
 \hat{\sigma}_{ij}=& \medspace \mathbbm{i} \mathrm{Kn} \big[ \omega_j(\hat{v}_i+\alpha_0\hat{q}_i)+\omega_i(\hat{v}_j+\alpha_0\hat{q}_j)\nonumber
\\
&-\omega_k(\hat{v}_k+\alpha_0\hat{q}_k)\delta_{ij}\big],
 \\
 \label{ccr2_ft}
 \hat{q}_i =& \medspace \mathbbm{i}\frac{c_p \mathrm{Kn}}{\mathrm{Pr}}\left(\omega_i \hat{T}+\alpha_0 \omega_j\hat{ \sigma}_{ij}\right).
\end{align}
Using Eqs.~\eqref{mass_ft} and \eqref{energy_ft},  Eq.~\eqref{ccr1_ft} simplifies to
\begin{align}
\label{ss1}
 \hat{\sigma}_{ij}=& \medspace\mathbbm{i} \mathrm{Kn}\big[ \omega_j(\hat{v}_i+\alpha_0\hat{q}_i)+\omega_i(\hat{v}_j+\alpha_0\hat{q}_j)\big]
\nonumber\\
 &+\mathrm{Kn}(h+\alpha_0g)\delta_{ij}.
\end{align}
Multiplying the above equation with  $\omega_j$ and $\omega_i \omega_j$,  we obtain
\begin{align}
 \label{ss11}
 \omega_j\hat{\sigma}_{ij}&=\mathbbm{i} \mathrm{Kn}[\omega^2(\hat{v}_i+\alpha_0\hat{q}_i)],
 \\
 \label{ss12}
 \omega_i \omega_j\hat{\sigma}_{ij}&=-\mathrm{Kn}\, \omega^2(h+\alpha_0 g),
\end{align}
respectively,  where $\omega_i\omega_i=|\omega_i|^2=\omega^2$.
Multiplying Eq.~\eqref{ccr2_ft} with $\omega_i$ and exploiting Eqs.~\eqref{energy_ft} and \eqref{ss12},  we obtain
\begin{align}
\label{theta_ft}
\hat{T}=\frac{g\,\mathrm{Pr}}{\omega^2 c_p \mathrm{Kn}}+\alpha_0 \mathrm{Kn}(h+\alpha_0 g).
\end{align}
Again,  multiplying Eq.~\eqref{mom_ft} with $\omega_i$ and exploiting Eq.~\eqref{ss12},  we obtain
\begin{align}
\label{p_ft}
\hat{p}&=\mathbbm{i}\frac{f_i\omega_i}{\omega^2}+\mathrm{Kn}(h+\alpha_0 g).
\end{align}
Now,  from Eqs.~\eqref{mom_ft} and \eqref{p_ft},  one can easily write
\begin{align}
\label{ss5}
\omega_j\hat{\sigma}_{ij}&=\mathbbm{i} f_i-\mathbbm{i}\omega_i\frac{f_k\omega_k}{\omega^2}-\omega_i \mathrm{Kn}(h+\alpha_0 g).
\end{align}
Substituting the value of $\hat{T}$ from Eq.~\eqref{theta_ft} and the value of $\omega_j\hat{\sigma}_{ij}$ from Eq.~\eqref{ss5} into Eq.~\eqref{ccr2_ft},  we obtain
\begin{align}
\label{q_ft}
\hat{q}_i&=\mathbbm{i} g\frac{\omega_i}{\omega^2}-\alpha_0\frac{c_p \mathrm{Kn}}{\mathrm{Pr}}f_k\left(\delta_{ik}-\frac{\omega_i\omega_k}{\omega^2}\right).
\end{align}
Now,  from Eqs.~\eqref{ss11}, \eqref{ss5} and \eqref{q_ft},
\begin{align}
\label{v_ft}
\hat{v}_i=& \medspace \frac{f_k}{\mathrm{Kn}}\left(\frac{\delta_{ik}}{\omega^2}-\frac{\omega_i\omega_k}{\omega^4}\right)
\nonumber\\
&+\alpha_0^2\frac{c_p \mathrm{Kn}}{\mathrm{Pr}}f_k\left(\delta_{ik}-\frac{\omega_i\omega_k}{\omega^2}\right)+\mathbbm{i} h\frac{\omega_i}{\omega^2}.
\end{align}
Finally,  using Eqs.~\eqref{v_ft} and \eqref{q_ft} in \eqref{ccr1_ft},  we obtain
\begin{align}
\label{sigma_ft}
\hat{\sigma}_{ij}=&\mathbbm{i} f_k
\left(\frac{\omega_j \delta_{ik}+\omega_i\delta{jk}}{\omega^2}-2\frac{\omega_i\omega_j\omega_k}{\omega^4}\right)\nonumber\\&-\mathrm{Kn}\left(\frac{\omega_i\omega_j}{\omega^2}-\frac{\delta_{ij}}{2}\right)(h+\alpha_0 g).
\end{align}
Applying the inverse Fourier transform in \eqref{theta_ft},  \eqref{p_ft} and \eqref{q_ft}--(\ref{sigma_ft}) with the help of formulae derived in Appendix\,\ref{app:A},  the field variables turn out to be
\begin{align}
\label{v}
v_i=&\frac{f_k}{\mathrm{Kn}}\left(\frac{x_ix_k}{4\pi r^2}-\frac{2\ln r-1}{8\pi}\delta_{ik}\right)\nonumber
\\&+\alpha_0^2\frac{c_p \mathrm{Kn}}{2\pi \mathrm{Pr}}f_k\left(\frac{2x_ix_k}{r^4}-\frac{\delta_{ik}}{r^2}\right)+\frac{h x_i}{2\pi r^2},
\\
\label{q}
q_i=&\frac{g x_i}{2\pi r^2}-\alpha_0\frac{c_p \mathrm{Kn}}{2\pi \mathrm{Pr}}f_k\left(\frac{2x_ix_k}{r^4}-\frac{\delta_{ik}}{r^2}\right),
\\
\label{p}
p=&\frac{f_i x_i}{2\pi r^2},
\\
\label{theta}
T=&-\frac{g\,\mathrm{Pr}\,\ln r}{2\pi\,\mathrm{Kn}\,c_p},
\\
\label{sigma}
\sigma_{ij}=&\frac{f_k x_k+2\mathrm{Kn}(h+\alpha_0 g)}{2\pi}\left(\frac{2x_ix_j}{r^4}-\frac{\delta_{ij}}{r^2}\right),
\end{align}
where $r=|x_i|$.
The field variables in \eqref{v}--\eqref{sigma} are the fundamental solutions of the linearized CCR model in 2D.
These fundamental solutions in the vectorial/tensorial notation are written as
\begin{align}
\label{vel}
\bm{v}(\bm{r})&=\frac{1}{8\pi \mathrm{Kn}}\bm{f}\cdot\bm{J}(\bm{r})+\frac{h\,\bm{r}}{2\pi r^2}+\frac{c_p \mathrm{Kn}}{2\pi \mathrm{Pr}}\alpha_0^2 \bm{f}\cdot \bm{K}(\bm{r}),\!
\\
\label{pressure}
p(\bm{r})&=\frac{\bm{f}\cdot\bm{r}}{2\pi r^2},
\\
\label{stress}
\bm{\sigma}(\bm{r})&=\frac{2 \mathrm{Kn}(h+ g\alpha_0)+\bm{f}\cdot\bm{r}}{2\pi} \bm{K}(\bm{r}),
\\
\label{temp}
T(\bm{r})&=-\frac{\mathrm{Pr}\, g}{2\pi \mathrm{Kn}\,c_p}\ln r,
\\
\label{heat}
\bm{q}(\bm{r})&=\frac{g}{2\pi}\frac{\bm{r}}{{r}^2}-\frac{c_p \mathrm{Kn}}{2\pi \mathrm{Pr}}\alpha_0 \bm{f}\cdot \bm{K}(\bm{r}),
\end{align}
where $r=|\bm{r}|$ and
\begin{align}
\bm{J}(\bm{r})&=\frac{2\bm{rr}}{r^2}-(2\ln r-1)\bm{I},
\\
\bm{K}(\bm{r})&=\frac{2\bm{rr}}{r^4}-\frac{\bm{I}}{r^2}.
\end{align}
It is worthwhile noticing that the fundamental solutions for the linearized NSF and G13 equations in 2D can be obtained directly from Eqs.~\eqref{vel}--\eqref{heat} by taking $\alpha_0 = 0$ and $\alpha_0 = 2/5$, respectively.
The fundamental solutions \eqref{vel}--\eqref{heat} need to be implemented with appropriate boundary conditions for a given problem.
We shall discuss their implementation and validation in Sec.\,\ref{sec:prob1} and Sec.\,\ref{sec:prob2} for two different problems.

\section{\label{sec:system}Boundary discretization}
We describe the construction of a system of algebraic equations for applying the MFS through the problem of flow past a complex geometry as depicted in Fig.~\ref{fig:arb}.
As an example, the geometry of the object in Fig.~\ref{fig:arb} is mathematically defined in the parametric form as
\begin{align}
\label{parametric}
(x,y)=\left(\frac{5}{4} a \cos{\theta}, \frac{1}{4}a (5 - \cos{5\theta}) \sin{\theta}\right)
\end{align}
with $0\leq \theta\leq 2\pi$ and $a\leq 1$ being the dilation factor.
\begin{figure}[t!]
\centering
\includegraphics[width=0.9\columnwidth]{figures/arb_schematic1.eps}
\caption{\label{fig:arb} Schematic diagram of a flow past an object of an arbitrary shape depicting the
boundary discretization and the placement of singularities outside the flow domain.
The red and blue arrows at each boundary node depict the normal (pointing toward the flow domain) and tangential directions at that node, respectively.}
\end{figure}
In the MFS, it is quite natural to place the singularity points outside the flow domain \cite{KA1964}.
Notwithstanding,
the location of the singularity points is a major concern as the results obtained from the MFS are highly sensitive toward the location of singularities \cite{Alves2009, CKL2016, CH2020}.
There are two most common ways of distributing singularities in the MFS.
One way is to place the singularities on a fictitious boundary of a very simple shape---irrespective of the shape of the object---with just one parameter to control; for example, on a circle in the two-dimensional case and on a sphere in the three-dimensional case, and the radius of the circle or sphere would be the controlling parameter.
Another way is to
recreate a dilated (or shrunk) fictitious boundary, which has the same shape as the boundary of the original object and to place the singularities on this fictitious boundary \cite{CH2020, LCK2013}---similarly to that shown in Fig.~\ref{fig:arb} as well.
The latter is also easy if the original boundary of the object can be described by a set of parametric equations having only a single controlling parameter, the dilation factor.
For the problem depicted in Fig.~\ref{fig:arb}, we have taken the fictitious boundary to be of the same shape as the original boundary.

Let $B$ be the number of the discretized boundary nodes and $S$ the number of singularity points. The boundary nodes and the singularities are placed at equispaced angles $\theta$ on the original and the fictitious boundary, respectively, and the distance between both boundaries can be varied by changing the value of the dilation factor $a$.
It may be noted that singularities need not be placed at equispaced angles in principle; nonetheless, we have done so for the sake of simplicity.
Let $\bm{x}_i^s$ and $\bm{x}_j^b$  be the position vectors of the $i^\mathrm{th}$ singularity site and the $j^\mathrm{th}$ boundary node, respectively,
then the position vector from the $i^\mathrm{th}$ singularity site to any position $\bm{x}$ in the domain is $\bm{r}_i=\bm{x}-\bm{x}_i^s$ and the position vector from the $i^\mathrm{th}$ singularity site to the $j^\mathrm{th}$ boundary node is $ \bm{r}_{ij}=\bm{x}_j^b-\bm{x}_i^s$.
It is important to note that the subscripts `$i$' and `$j$' are now being used for denoting the $i^\mathrm{th}$ singularity site and $j^\mathrm{th}$ boundary node and consequently, the repetition of indices henceforth shall \emph{not} imply the Einstein summation per se, unless stated otherwise (particularly, in Appendix~\ref{app:A}, wherein the Einstein summation does hold over the repeated indices).
Since the point sources $\bm{f}$, $g$ and $h$---of different strengths---are to be put at each singularity site, there are four degrees of freedom corresponding to each singularity point [two scalars $g$ and $h$ from the point heat and mass sources, and two components $f_1$ and $f_2$ of the point force vector $\bm{f}=(f_1,f_2)^\mathsf{T}$].
In total, we have $4\times S$ unknowns, which are determined typically by satisfying the boundary conditions at the boundary points.
Once the location of the singularity points is decided, the next step in the implementation of the MFS is superposition of the fundamental solutions associated with each singularity site, which makes sense because of the linearity of equations and gives the value of the field variables at the $j^\mathrm{th}$ boundary node.
Superimposing the fundamental solutions \eqref{vel}--\eqref{heat} for each singularity site, the field variables at the $j^\mathrm{th}$ boundary node read
\begingroup
\allowdisplaybreaks
\begin{align}
\label{vj}
\bm{v}_
&=\sum_{i=1}^{S}\bigg[\frac{\bm{f}_i\cdot\bm{J}(\bm{r}_{ij})}{8\pi \mathrm{Kn}}+\frac{h_i}{2\pi}\frac{\bm{r}_{ij}}{r_{ij}^2}+\frac{c_p \mathrm{Kn}}{2\pi \mathrm{Pr}}\alpha_0^2 \bm{f}_i\cdot \bm{K}(\bm{r}_{ij})\bigg],
\\
\label{Pj}
p_
&=\sum_{i=1}^{S}\frac{\bm{f}_i\cdot \bm{r}_{ij}}{2\pi r_{ij}^2},
\\
\bm{\sigma}_
&=\sum_{i=1}^{S}\frac{2 \mathrm{Kn}\,(h_i+g_i\,\alpha_0)+\bm{f}_i\cdot \bm{r}_{ij}}{2\pi} \bm{K}(\bm{r}_{ij}),
\\
\label{Tj}
T_
&=-\sum_{i=1}^{S}\frac{ g_i\, \mathrm{Pr}}{c_p \mathrm{Kn}}\frac{\ln r_{ij}}{2\pi},
\\
\label{qj}
\bm{q}_
&=\sum_{i=1}^{S} \bigg[\frac{g_i}{2\pi}\frac{\bm{r}_{ij}}{r_{ij}^2}-\frac{c_p \mathrm{Kn}}{2\pi \mathrm{Pr}}\alpha_0 \bm{f}_i\cdot \bm{K}(\bm{r}_{ij}) \bigg],
\end{align}
\endgroup
where $r_{ij}=|\bm{r}_{ij}|$; $\bm{f}_i=({f_1}_i,{f_2}_i)^\mathsf{T}$, $g_i$ and $h_i$ are the point force (vector), point heat source and point mass source, respectively, applied on the $i^\mathrm{th}$ singularity site; and
\begin{align}
\bm{J}(\bm{r}_{ij})=&\frac{2\bm{r}_{ij}\bm{r}_{ij}}{r_{ij}^2}-(2\ln r_{ij}-1)\bm{I},
\\
\bm{K}(\bm{r}_{ij})=&\frac{2\bm{r}_{ij}\bm{r}_{ij}}{r_{ij}^4}-\frac{\bm{I}}{r_{ij}^2}.
\end{align}
This system is solved for the unknowns ${f_1}_i,{f_2}_i,g_i,h_i$, $i\in\{1,2,3,\dots,S\}$ by employing the boundary conditions at each boundary node.
Once the unknowns ${f_1}_i,{f_2}_i,g_i,h_i$ for $i\in\{1,2,3,\dots,S\}$ are found,
the flow variables at any position $\bm{x}$ in the flow domain can be determined
simply by dropping the subscript `$j$' everywhere
in Eqs.~\eqref{vj}--\eqref{qj}.
For instance, the velocity $\bm{v}\equiv\bm{v}(\bm{x})$ at a position $\bm{x}$ in the flow domain is given by
\begin{align}
\bm{v}&=\sum_{i=1}^{S}\bigg[\frac{\bm{f}_i\cdot\bm{J}(\bm{r}_{i})}{8\pi \mathrm{Kn}}+\frac{h_i}{2\pi}\frac{\bm{r}_{i}}{r_{i}^2}+\frac{c_p \mathrm{Kn}\alpha_0^2}{2\pi \mathrm{Pr}} \bm{f}_i\cdot \bm{K}(\bm{r}_{i})\bigg].\!
\end{align}
The other flow variables are obtained from Eqs.~\eqref{Pj}--\eqref{qj} analogously.
The above procedure to evaluate flow variables
works for any geometry and we have implemented this in a numerical framework.
We shall elaborate on the placement of boundary nodes and source points, formation and solution of the system separately corresponding to the two problems in Sec.\,\ref{sec:prob1} and Sec.\,\ref{sec:prob2}.
\section{\label{sec:prob1}Vapor flow confined between two coaxial cylinders}
For the validation of the developed numerical framework, we revisit the problem of a rarefied vapor flow confined between two concentric cylinders.
The same problem was investigated by Onishi \cite{Onishi1977} with the linearized BGK model and the diffuse reflection boundary conditions.
\subsection{\label{prob1_disc}Problem description}
Let us consider a moderately rarefied vapor confined between the condensed phases of two concentric infinitely long circular cylinders of radii $\tilde{R}_1$ and $\tilde{R}_2$, where $\tilde{R}_1<\tilde{R}_2$.
Owing to the axial symmetry, it is sufficient to investigate the problem in 2D.
A cross-sectional (two-dimensional) view of the problem is illustrated in Fig.~\ref{fig:coaxial_schematic}.
\begin{figure}[!b]
\centering
\includegraphics[width=\columnwidth]{figures/coaxial.eps}
\caption{\label{fig:coaxial_schematic}Cross-sectional view of the vapor flow confined between two coaxial cylinders.}
\end{figure}
For the purpose of non-dimensionalization, we take the inner radius as the characteristic length $\tilde{L}$, i.e.~$\tilde{L}=\tilde{R}_1$.
Consequently, the dimensionless radii of the inner and outer cylinders are $r_1=\tilde{R}_1/\tilde{L}=1$ and $r_2=\tilde{R}_2/\tilde{L}$, respectively.
The condensed phases of the vapor at the inner and outer cylinders are assumed to be negligibly thin.
Let the temperatures of the inner and outer condensed phases be maintained at uniform temperatures $\tilde{T}_0$ and $\tilde{T}_s$, respectively.
Moreover, let the saturation pressures of the condensed phases corresponding to the temperatures $\tilde{T}_0$ and $\tilde{T}_s$ be $\tilde{P}_0$ and $\tilde{P}_s$, respectively; see Fig.~\ref{fig:coaxial_schematic}.
Again, for the purpose of linearization and non-dimensionalization, we take the temperature at the inner wall $\tilde{T}_0$  as the reference temperature and the saturation pressure at the inner wall $\tilde{P}_0$ as the reference pressure.
Thus, the dimensionless perturbations in the temperature and saturation pressure at the inner wall vanish, and the dimensionless perturbations in the temperature and saturation pressure at the outer wall read
\begin{align}
\tau_s=\frac{\tilde{T}_s-\tilde{T}_0}{\tilde{T}_0}
\quad\text{and}\quad
p_s=\frac{\tilde{P}_s-\tilde{P}_0}{\tilde{P}_0},
\end{align}
respectively.
\subsection{Analytic solution of Onishi \cite{Onishi1977}}
Onishi \cite{Onishi1977} investigated the problem by employing an asymptotic theory \cite{Sone2002}.

According to this theory,
a field variable $\tilde{h}$ of the gas can be written as
\begin{align}
\tilde{h}=\tilde{h}_H+\tilde{h}_K,
\end{align}
where $\tilde{h}_H$ is referred to as the hydrodynamic part or the Hilbert part that describes the flow behavior in the bulk of the domain and $\tilde{h}_K$ is referred to as the kinetic boundary layer part or the Knudsen layer part that can be seen as a correction to the Hilbert part and is significant only in a small layer near an interface.
Both $\tilde{h}_H$ and $\tilde{h}_K$ for all field variables are expanded in power series of the Knudsen number, and the contribution at each power of the Knudsen number is then computed by means of the considered model (the BGK model in [\onlinecite{Onishi1977}]) and appropriate boundary conditions (the diffuse reflection boundary conditions in [\onlinecite{Onishi1977}]).

The linearized CCR model is anyway not able to predict Knudsen layers. Therefore, it makes sense to compare the results obtained from the MFS only with the Hilbert part of the solution given in Ref.~[\onlinecite{Onishi1977}].
For the problem under consideration and for the linearized BGK model with the diffuse reflection boundary conditions, the Hilbert part of the solution is indeed straightforward to determine by solving a set of simple ordinary differential equations analytically, see Ref.~[\onlinecite{Onishi1977}].
Denoting the radius ratio by $\beta=r_2/r_1$ and the ratio of $p_s$ to $\tau_s$ by
$\gamma=p_s/\tau_s$, the analytic solution obtained from the linearized BGK model with the diffuse reflection boundary conditions for $\mathrm{Kn}\approx 0$ is given by \cite{Onishi1977}
\begin{align}
\label{analy_p}
p=& \medspace
p_s\left(\frac{1}{r_1}+\frac{1}{r_2}\right)^{-1}\frac{1}{r_1},
\\
\label{analy_vr}
v_r=&-\frac{p_s}{C_0}\left(\frac{1}{r_1}+\frac{1}{r_2}\right)^{-1}\frac{1}{r},
\\
\label{analy_T}
T=& \medspace
\tau_s\left[\left(1-\frac{D_0}{C_0}\gamma\right)\frac{\ln r}{\ln \beta}-\left(1-\frac{D_0}{C_0}\gamma\right)\frac{\ln r_1}{\ln \beta}\right]
\nonumber\\
&+\frac{D_0}{C_0}\gamma\tau_s\left(\frac{1}{r_1}+\frac{1}{r_2}\right)^{-1}\frac{1}{r_1},
\\
\label{analy_q}
q_r=&0,
\end{align}
where $C_0=2.132039$ and $D_0=0.4467494$.

\subsection{\label{prob1_bcs}Boundary conditions and implementation of the MFS}
We shall revisit the problem described above by means of the MFS applied on the linearized CCR model.
Recall that we have already determined the fundamental solutions of the linearized CCR model and outlined the way to implement them in Sec.\,\ref{sec:solutions} for a general two-dimensional object.
The solution for the field variables at the $j^\mathrm{th}$ boundary node  (Eqs.~\eqref{vel}--\eqref{heat}) can directly be used once the boundary nodes and singularity points for the present problem have been decided.

\begin{figure}[!b]
\centering
\includegraphics[width=\columnwidth]{figures/sing_coaxial.eps}
\caption{\label{fig:coaxial_discr}Schematic of the boundary nodes on the boundaries and singularity points outside the flow domain for the problem illustrated in Fig.~\ref{fig:coaxial_schematic}. The red and blue arrows at each boundary node depict the normal (pointing toward the flow domain) and tangential directions at that node, respectively.%
}
\end{figure}
Since the singularity sites are to be placed outside of the computational domain, we assume the source points to be placed on two fictitious circular boundaries, one inside the circle associated with the inner cylinder and the other outside the circle associated with the outer cylinder, as shown in Fig.~\ref{fig:coaxial_discr}.
Note that both fictitious boundaries are concentric with the circles associated with the cylinders.
Let the radii of the inner and outer fictitious boundaries be $\tilde{S}_1$ and $\tilde{S}_2$, respectively. %
For simplicity, we consider $N_s$ equispaced source points on each of the two fictitious boundaries and $N_b$ equispaced boundary nodes on each of the actual boundaries (the boundaries of the inner and outer cylinders).
As explained in Sec.\,\ref{sec:system}, we have four degrees of freedom corresponding to each source point, and the total number of singularity points for the problem under consideration is $S=2N_s$. Thus, there will be a total of $4\times S = 4\times2N_s$ unknowns in the problem.
Accordingly, the summations in Eqs.~\eqref{vel}--\eqref{heat} will run from $i=1$ to $i=2N_s$.

Boundary conditions at the $j^{\mathrm{th}}$ boundary node are obtained from \eqref{bc_1}--\eqref{bc_3} by replacing the flow variables and the normal and tangent vectors with their respective values at the $j^{\mathrm{th}}$ boundary node.
Furthermore, since the walls of the cylinders are fixed, $\bm{v}^I=\bm{0}$.
Consequently, the boundary conditions at the $j^{\mathrm{th}}$ boundary node read
\begin{align}
\label{bc_11}
\bm{v}_j\cdot \bm{n}_j=& -\eta_{11}(p_j-p_{\mathrm{sat}} + \bm{n}_j\cdot\bm{\sigma}_j\cdot\bm{n}_j)\nonumber
\\&+\eta_{12}(T_j-T^I+\alpha_0 \bm{n}_j\cdot\bm{\sigma}_j\cdot\bm{n}_j),
\\[2ex]
\label{bc_21}
\bm{q}_j\cdot \bm{n}_j=& \medspace
\eta_{12}(p_j-p_{\mathrm{sat}} + \bm{n}_j\cdot\bm{\sigma}_j\cdot\bm{n}_j)\nonumber\\&-(\eta_{22}+2\tau_0)(T_j-T^I+\alpha_0 \bm{n}_j\cdot\bm{\sigma}_j\cdot\bm{n}_j),
\end{align}
\begin{align}
\label{bc_31}
\bm{t}_j\cdot\bm{\sigma}_j\cdot\bm{n}_j=&-\varsigma (\bm{v}_j+\alpha_0 \bm{q}_j)\cdot\bm{t}_j.
\end{align}
The dimensionless perturbations in saturation pressures at the inner and outer interfaces are $p_{\mathrm{sat}}=0$ and $p_{\mathrm{sat}}=p_s$, respectively, and the dimensionless perturbations in temperatures at the inner and outer interfaces are $T^I=0$ and $T^I=\tau_s$, respectively, which need to be replaced in boundary conditions \eqref{bc_11}--\eqref{bc_31} accordingly.
Note that boundary conditions \eqref{bc_11}--\eqref{bc_31} are to be satisfied at $B=2N_b$ boundary nodes.
On substituting the values of the field variables at the $j^\mathrm{th}$ boundary node from \eqref{vj}--\eqref{qj} into boundary conditions \eqref{bc_11}--\eqref{bc_31}, the resulting system of equations (associated with the $j^\mathrm{th}$ boundary node) can be written in a matrix form as
\begin{align}
\label{matrix}
\sum_{i=1}^{S} M_{ji} \bm{u}_i=\bm{b}_j,
\end{align}
for the unknown vector associated with the $i^\mathrm{th}$ singularity $\bm{u}_i=({f_1}_i,{f_2}_i,g_i,h_i)^\mathsf{T}$.
Here, $M_{ji}$'s are coefficient matrices of dimensions  $3\times 4$ and $\bm{b}_j$ is the $3\times 1$ vector containing the interface properties, such as $p_s$ and $\tau_s$.
We collect all such $B$ systems into a new system
\begin{align}
\label{collocation_system}
\mathcal{M} \bm{\mathcal{X}}=\bm{\mathcal{B}},
\end{align}
where $\bm{\mathcal{X}} = \big({f_1}_1,{f_2}_1,g_1,h_1,{f_1}_2,{f_2}_2,g
_2,h_2,\dots,{f_1}_S,{f_2}_S,\allowbreak g_S,h_S\big)^\mathsf{T}$ is the vector containing all $4S$ unknowns, the matrix $\mathcal{M}$---containing the coefficients of the unknowns---has dimensions $3B\times 4S$ (or $6N_b\times 8N_s$) and is referred to as the collocation matrix.
We have solved system \eqref{collocation_system} in the computer algebra software, Mathematica\textsuperscript{\textregistered} using the method of least squares.
For the identification purpose, the first $N_s$ singularity points ($i=1,2,\dots,N_s$) in our code belong to the inner fictitious boundary and the rest $N_s$ singularity points ($i=N_s+1,N_s+2,\dots,2N_s$) to the outer fictitious boundary.
Similarly, the first $N_b$ boundary nodes ($j=1,2,\dots,N_b$) belong to the actual inner boundary and the rest $N_b$ boundary nodes ($j=N_b+1,N_b+2,\dots,2N_b$) to the actual outer boundary.
\subsection{Results and discussion}
For numerical computations, we have taken $N_b=100$ boundary nodes on each of the actual boundaries and $N_s=100$ singularity points on each of the fictitious boundaries.
The dimensionless radii of the original and fictitious boundaries are taken as $r_1=1$, $r_2=2$, $s_1=\tilde{S}_1/\tilde{R}_1=0.5$ and $s_2=\tilde{S}_2/\tilde{R}_2=4$.

Figure~\ref{fig:temp_gamma} illustrates the variation of the (scaled) temperature of the vapor in the radial direction for $\mathrm{Kn} \approx 0$ and for different values of the parameter $\gamma \,(=p_s/\tau_s)$, wherein $\tau_s=4$ is fixed and $p_s$ is being varied for varying $\gamma$.
The solid lines represent the results obtained from our numerical framework based on the MFS while the symbols delineate the results from Eq.~\eqref{analy_T}, which was obtained analytically for $\mathrm{Kn} \approx 0$ through an asymptotic theory \cite{Sone2002} performed on the linearized BGK model in Ref.~[\onlinecite{Onishi1977}].
It is evident from the figure that the results obtained with the MFS in the present work are in an excellent agreement with the analytic results  from the linearized BGK model for $\mathrm{Kn} \approx 0$.\begin{figure}[!b]
\centering
\includegraphics[width=\columnwidth]{figures/k0_validation}
\caption{\label{fig:temp_gamma}
Variation of the (scaled) temperature in the gap between the two cylinders for different values of $\gamma$.
The solid lines denote the results obtained from the MFS applied on the CCR model and the symbols indicate the analytic results from Eq.~\eqref{analy_T}, which was obtained analytically from the linearized BGK model for $\mathrm{Kn} \approx 0$ in Ref.~[\onlinecite{Onishi1977}].
The other parameters are $N_b=100$, $N_s=100$, $r_1=1$, $r_2=2$, $s_1=0.5$, $s_2=4$.
}
\end{figure}
Although not shown here for brevity, the results for the pressure and velocity from the MFS are also in excellent agreement with the analytic results from Eqs.~\eqref{analy_p} and \eqref{analy_vr} for $\mathrm{Kn} \approx 0$.
It is also evident from Fig.~\ref{fig:temp_gamma} that the temperature increases on moving away from the inner cylinder toward the outer cylinder for smaller values of $\gamma$ (red and blue curves and symbols in the figure) and vice versa for larger values of $\gamma$ (green and magenta curves and symbols in the figure).
This indicates the existence of a reverse temperature gradient after a  critical value of $\gamma$.
Indeed, at this critical value of $\gamma$, the (scaled) temperature remains constant along the radial direction.
An expression for this critical value of $\gamma$ from the asymptotic theory \cite{Sone2002} is given by \cite{Onishi1977}
\begin{align}
\label{gammac}
\gamma_c=\frac{C_0}{D_0}\left[1-\mathrm{Kn} \frac{C_0}{D_0}(0.124226)\left(\frac{1}{r_1}-\frac{1}{r_2}\right) + \mathcal{O}(\mathrm{Kn}^2)\right].
\end{align}
For $\mathrm{Kn} \approx 0$, the critical value of $\gamma$ from the above expression is $\gamma_c = C_0/D_0 \approx 4.772337$.
From the MFS presented here, the critical value of $\gamma$ for $\mathrm{Kn} \approx 0$ turns out to be $\gamma_c \approx 4.7723$, which is also very close to that computed from the above expression.
The phenomenon of reverse temperature gradient can be understood form boundary condition \eqref{bc_21} as follows.
There are two factors determining the normal heat flux component in boundary condition \eqref{bc_21} according to which the evaporation/condensation rate depends on (i) the difference between the pressure and saturation pressure, and (ii) the temperature difference between the temperatures of the gas (or vapor) and and interface.
The temperature gradient gets reversed when one dominates the other.

To examine the capabilities of the developed method, we also study the problem for higher Knudsen numbers.
Figure~\ref{fig:temp_gamma_Kn0p1} exhibits the variation of the (scaled) temperature of the vapor in the radial direction for $\mathrm{Kn} = 0.1$ and for different values of the parameter $\gamma$.
The  solid lines again represent the results obtained from our numerical framework based on the MFS but the symbols now denote the data from the linearized BGK model taken directly from Ref.~[\onlinecite{Onishi1977}].
\begin{figure}[!b]
\centering
\includegraphics[width=\columnwidth]{figures/k01_validation}
\caption{\label{fig:temp_gamma_Kn0p1}Same as Fig.~\ref{fig:temp_gamma} except for the symbols denote the data from Ref.~[\onlinecite{Onishi1977}] for $\mathrm{Kn} = 0.1$ obtained using the linearized BGK model.
}
\end{figure}
It is clear from the figure that the results from the MFS are in good agreement with those from the  linearized BGK model even for $\mathrm{Kn} = 0.1$; nonetheless, the quantitative differences in the results from both methods are now noticeable.
In addition, figure~\ref{fig:temp_gamma_Kn0p1} also shows the existence of a reverse temperature gradient.
For $\mathrm{Kn} = 0.1$, the critical value of $\gamma$, at which the phenomenon of reverse temperature gradient occurs, computed from the MFS
is $\gamma_c=4.66247$ whereas its reported value from the linearized BGK model in Ref.~[\onlinecite{Onishi1977}] is $\gamma_c=4.63087$.

\begin{figure}[!b]
\centering
\includegraphics[width=\columnwidth]{figures/heat_coaxial.eps}
\caption{
\label{fig:heat_coaxial}
Variation of the (scaled) radial heat flux with $\gamma$.
The solid lines denote the results obtained from the MFS applied on the CCR model and the symbols indicate the data taken directly from Ref.~[\onlinecite{Onishi1977}], which were obtained using the linearized BGK model.
The other parameters are the same as those for Fig.~\ref{fig:temp_gamma}.
}
\end{figure}

To have further insight on the reverse temperature gradient, the (scaled) radial heat flux at the actual inner boundary (i.e.~at $r=1$) is plotted against $\gamma$ in Fig.~\ref{fig:heat_coaxial}.
The solid lines and symbols denote the results from the MFS in the present work and the data from the linearized BGK model given in Ref.~[\onlinecite{Onishi1977}], respectively.
It is apparent from the figure that our results for the radial heat flux are also in good agreement with the data from the linearized BGK model for a smaller value of the Knudsen number ($\mathrm{Kn}=0.1$ in the figure); however, for a higher value of the Knudsen number ($\mathrm{Kn}=0.2$ in the figure), there is a noticeable mismatch between the results obtained from the MFS and the data from the linearized BGK model given in Ref.~[\onlinecite{Onishi1977}].
A plausible reason for this discrepancy could be the truncation of power series at the first order  in  Ref.~[\onlinecite{Onishi1977}] because the neglected terms in the series could have significant contributions for larger values of the Knudsen number.
Figure~\ref{fig:heat_coaxial} also shows that for each value of the Knudsen number, there is a $\gamma$ at which the radial heat flux changes its sign.
This $\gamma$ is indeed the same as the $\gamma_c$ described above, at which reversal of the temperature gradient takes place.

Through the plots of heat flux lines, also not shown here for brevity, it has been found that, in the case of $\tau_s>0$, heat  flows from the outer cylinder toward the inner cylinder for $\gamma<\gamma_c$ and vice versa for $\gamma>\gamma_c$ .
This makes sense in view of Figs.~\ref{fig:temp_gamma} and \ref{fig:temp_gamma_Kn0p1}.
The direction of heat flow reverses in both cases when $\tau_s$ is taken to be negative or, in other words, when the initial temperature of the inner cylinder is taken higher than that of the outer cylinder.

Figure~\ref{fig:vel_coaxial} displays the (scaled) radial velocity plotted against $\gamma$ for $\mathrm{Kn} \approx 0$ and $\mathrm{Kn} =0.1$.
The solid lines are again the results from the MFS in the present work while the symbols in the case of $\mathrm{Kn} \approx 0$ denote the results from Eq.~\eqref{analy_vr} and those in the case of $\mathrm{Kn} = 0.1$ denote the data taken from Ref.~[\onlinecite{Onishi1977}]; nevertheless, in both cases symbols denote the results from the linearized BGK model.
\begin{figure}[!t]
\includegraphics[width=\columnwidth]{figures/vel_coaxial.eps}
\caption{\label{fig:vel_coaxial}Variation of the (scaled) radial velocity with $\gamma$.
The solid lines denote the results obtained from the MFS applied on the CCR model and the symbols indicate those from the linearized BGK model (from Eq.~\eqref{analy_vr} in the case of $\mathrm{Kn} \approx 0$ and directly from Ref.~[\onlinecite{Onishi1977}] in the case of  $\mathrm{Kn} = 0.1$).
The other parameters are the same as those for Fig.~\ref{fig:temp_gamma}.}
\end{figure}
The figure also demonstrates a good agreement between the results from the method developed in the present work and those from the linearized BGK model.

\section{\label{sec:prob2}Rarefied gas flow between two non-coaxial cylinders}
In this section, we revisit the problem of flow induced by a temperature difference in a rarefied gas confined between two non-coaxial cylinders via the MFS developed above.
The same problem was investigated numerically by Aoki, Sone and Yano \cite{Aoki_Sone_1989} with the linearized BGK model and the diffuse reflection boundary conditions.
\subsection{Problem description}
Let us consider a rarefied (monatomic) gas confined between two infinitely long circular cylinders of radii $\tilde{R}_1$ and $\tilde{R}_2$ (with $\tilde{R}_1<\tilde{R}_2$) that are not coaxial.
Again, owing to the axial symmetry, it is sufficient to investigate the problem in 2D.
\begin{figure}[!b]
\centering
\includegraphics[width=\columnwidth]{figures/noncoaxial} \caption{\label{fig:noncoaxial_schematic}Cross-sectional view of the flow of a rarefied gas confined between two non-coaxial cylinders having different wall temperatures.}
\end{figure}
Let the locations of both cylinders be fixed according to the cross-sectional view portrayed in Fig.~\ref{fig:noncoaxial_schematic} so that the centers of the circles associated with the inner and outer cylinders be at the origin and at $(0, - \tilde{d})$, respectively.
Furthermore, let the temperatures of the inner and outer cylinders be kept fixed at $\tilde{T}_i = \tilde{T}_0$ and $\tilde{T}_o = \tilde{T}_0 (1+\Delta\tau)$, respectively, with $\Delta\tau$ being sufficiently small in comparison to $\tilde{T}_0$ so that the linear theory remains meaningful.

For the purpose of non-dimensionalization, we again take the radius of the inner cylinder as the characteristic length $\tilde{L}$, i.e.~$\tilde{L}=\tilde{R}_1$.
Consequently, the dimensionless radii of the inner and outer cylinders are $r_1=\tilde{R}_1/\tilde{L}=1$ and $r_2=\tilde{R}_2/\tilde{L}$, respectively, and the dimensionless distance between the centers of the cylinders is $d =\tilde{d}/\tilde{L}$.
Furthermore, for the purpose of the linearization and non-dimensionalization, the equilibrium pressure of the gas $\tilde{p}_0$ is taken as the reference pressure and the temperature of the inner cylinder $\tilde{T}_i$ as the reference temperature so that the dimensionless perturbations in temperatures on the inner and outer walls are $T_i = (\tilde{T}_i - \tilde{T}_i) / \tilde{T}_i=0$ and $T_o = (\tilde{T}_o - \tilde{T}_i) / \tilde{T}_i = \Delta \tau$, respectively.
For comparing the results from the present method with those of Ref.~[\onlinecite{Aoki_Sone_1989}], the parameters are fixed to $r_2 = 2$, $d = 0.5$ and $\Delta \tau = 1$.

\begin{figure}[!b]
\centering
\includegraphics[width=\columnwidth]{figures/sing_noncoaxial.eps}
\caption{\label{fig:noncoaxial_discr}Schematic of the boundary nodes on the boundaries and singularity points outside the flow domain for the problem illustrated in Fig.~\ref{fig:noncoaxial_schematic}.
The red and blue arrows at each boundary node depict the normal (pointing toward the flow domain) and tangential directions at that node, respectively.%
}
\end{figure}
\subsection{Boundary Conditions and implementation of the MFS}
In order to place the singularity sites outside the computational domain, we again assume the source points to be placed on two fictitious circular boundaries, one inside the circle associated with the inner cylinder and the other outside the circle associated with the outer cylinder, as shown in Fig.~\ref{fig:noncoaxial_discr}.
The inner (outer) fictitious boundary is concentric with the circle associated with the inner (outer) cylinder.
Let the radii of the inner and outer fictitious boundaries be $\tilde{S}_1$ and $\tilde{S}_2$, respectively.
Consequently, the dimensionless radii of the inner and outer fictitious boundaries are $s_1 = \tilde{S}_1/\tilde{L}$ and $s_2 = \tilde{S}_2/\tilde{L}$.
Similarly to the above, we consider $N_s$ equispaced source points on each of the two fictitious boundaries and $N_b$ equispaced boundary nodes on each of the actual boundaries (the boundaries of the inner and outer cylinders).

Since the walls of the cylinders are fixed for this problem as well, $\bm{v}^I=\bm{0}$.
Hence, the boundary conditions \eqref{bc_11}--\eqref{bc_31} at the $j^{\mathrm{th}}$ boundary node hold true for the present problem as well.
However, since the present problem does not involve evaporation and condensation, the evaporation/condensation coefficient $\vartheta$ is zero for this problem.
Consequently, boundary conditions \eqref{bc_11}--\eqref{bc_31} for the problem under consideration further reduce to
\begin{align}
\label{bc1}
\bm{v}_j\cdot \bm{n}_j&= 0,
\\
\label{bc2}
\bm{q}_j\cdot \bm{n}_j&= -2\tau_0(T_j-T^I+\alpha_0 \,\bm{n}_j\cdot\bm{\sigma}_j\cdot\bm{n}_j),
\\
\label{bc3}
\bm{t}_j\cdot\bm{\sigma}_j\cdot\bm{n}_j&=-\varsigma (\bm{v}_j +\alpha_0^{\prime} \bm{q}_j)\cdot\bm{t}_j.
\end{align}
Note that the coefficient $\alpha_0$ in boundary condition \eqref{bc3} has been changed to $\alpha_0^{\prime}=1/5$ (see, e.g.~Refs.~[\onlinecite{Struchtrup2005, ST2008, TTS2009}]) in order to have a fair comparison with the findings of Ref.~[\onlinecite{Aoki_Sone_1989}].
The interface temperature $T^I$ in boundary condition \eqref{bc2} is $0$ for the inner cylinder and $\Delta \tau$ for the outer cylinder.

The construction of the collocation matrix and the formation of system \eqref{collocation_system} for the present problem is exactly similar to that demonstrated in Sec.~\ref{prob1_bcs}.
We have again solved system  \eqref{collocation_system} for the present problem analogously in the computer algebra software, Mathematica\textsuperscript{\textregistered} using the method of least squares to determine the unknowns ${f_1}_1,{f_2}_1,g_1,h_1,{f_1}_2,{f_2}_2,g
_2,h_2,\dots,{f_1}_S,{f_2}_S,\allowbreak g_S,h_S$.
\subsection{Results and discussion}
We have computed the results numerically by taking the parameters as $\Delta\tau=1$, $N_b = N_s = 100$, $r_1=1$, $r_2=2$, $s_1=0.5$, $s_2=4$.
In all the figures below, the solid lines represent the results obtained with the MFS in the present work and the symbols denote the data taken from Ref.~[\onlinecite{Aoki_Sone_1989}], which were obtained using the linearized BGK model.

Figure~\ref{fig:vel} illustrates the variation of the tangential component of the (dimensionless) velocity on the right halves of the inner (top row) and outer (bottom row) circles associated with the respective cylinders with respect to the angle $\theta$, which is the angle measured from the negative $y$-axis anticlockwise around the center of the inner circle as shown in Fig~\ref{fig:noncoaxial_discr}.
The angle has been taken in this way in order to maintain the geometrical similarity with Ref.~[\onlinecite{Aoki_Sone_1989}].
The unit tangential directions on the inner and outer circles are marked in Fig~\ref{fig:noncoaxial_discr} with blue arrows.
Figure~\ref{fig:vel} shows that the tangential components of the velocity for both inner and outer circles remain zero at $\theta = 0$ and $\theta = \pi$ and that they attain the maximum values somewhere in $(0, \pi/2)$.
Furthermore, the value of $\theta$ at which the maximum is attained also shifts more toward $\theta = \pi/2$ on increasing the value of the Knudsen number.
Figure~\ref{fig:vel} evinces that the results from the MFS applied on the CCR model (solid lines) are in reasonably good agreement with those from the linearized BGK model for small Knudsen numbers (green lines and symbols) and that the differences between the results from both methods become more and more prominent with increasing Knudsen numbers (red and blue lines and symbols), where the present method starts overpredicting the results, though the trends from both methods remain qualitatively similar to each other even for high Knudsen numbers.
\begin{figure}[!tb]
\centering
\includegraphics[width=\columnwidth]{figures/inner_cyl}
\\[3ex]
\includegraphics[width=\columnwidth]{figures/outer_cyl}
\caption{\label{fig:vel}
Tangential velocity on the right halves of the inner and outer circles associated with the respective cylinders plotted against the angle $\theta$ for different values of the Knudsen number and for $\Delta \tau = 1$.
The solid lines denote the results obtained from the MFS applied on the CCR model and the symbols indicate the data from the linearized BGK model \cite{Aoki_Sone_1989}.
The other parameters are the same as those for Fig.~\ref{fig:temp_gamma}.
}
\end{figure}

Figure~\ref{fig:vel}, in other words,  also reveals that at $\theta = 0$ and $\theta = \pi$ the flow can happen only in the normal directions.
This prompts us to draw the flow streamlines in Fig.~\ref{fig:stream}.
\begin{figure}[htb]
\centering
\includegraphics[width=\columnwidth]{figures/stream.eps}
\caption{\label{fig:stream}
Velocity streamlines in a cross section for the problem described in Sec.~\ref{sec:prob2} at $\mathrm{Kn}=0.1$.
The other parameters are the same as those for Fig.~\ref{fig:vel}.
}
\end{figure}
The streamlines in Fig.~\ref{fig:stream} shows that at the narrowest gap (near $\theta = 0$), the gas starts moving from the outer (hotter) cylinder toward the inner (colder) cylinder and flows along the surface of the inner cylinder on both halves until it reaches $\theta = \pi$, at which it can flow only in the normal direction.
Therefore, at the widest gap (near $\theta = \pi$), the gas flows from the inner cylinder toward the outer cylinder and returns back from there toward the narrowest gap along the surface of the outer cylinder (but in the opposite directions due to symmetry along the $y$-axis).
This renders two counter-directional circulating flows, one in the left half of the domain and the other in the right half of the domain.
The directions of the circulating flows reverse on taking $\Delta\tau<0$, i.e.~when the inner cylinder is at a higher temperature than the outer one.
From the considered values of the Knudsen number, the directions of the circulating flows apparently do not depend on the Knudsen number.

The superposition of the second components of the point force vectors at the inner source points gives the total drag force $D$ on the inner cylinder, i.e.%
\begin{align}
D=-\sum_{i=1}^{N_s} {f_2}_i,
\end{align}
where $ i=1,2,\dots,N_s$ refer to the points on the inner fictitious boundary and the negative sign represents the direction opposite  to the flow.
Variation of the drag force with the Knudsen number is illustrated in Figure~\ref{drag}, which shows good agreement between the results from the MFS applied on the CCR model (solid lines) and those from the linearized BGK model (symbols) even for high Knudsen numbers (especially, for $\mathrm{Kn}\lesssim 2$).
\begin{figure}[!tb]
\centering
\includegraphics[width=\columnwidth]{figures/drag.eps}
\caption{\label{drag}
Drag force on the inner cylinder plotted against the Knudsen number.
The solid lines denote the results obtained from the MFS applied on the CCR model and the symbols indicate the data from the linearized BGK model \cite{Aoki_Sone_1989}.
The other parameters are the same as those for Fig.~\ref{fig:vel}.
}
\end{figure}
This was actually not the case for tangential velocity displayed in Fig.~\ref{fig:vel}, where the differences between the results from the two models were noticeable for high Knudsen numbers.
This shows that the CCR model is capable of predicting the global quantities, e.g.~the drag force, quite accurately but is incapacitated of predicting the local quantities, e.g.~the velocity and temperature, for high Knudsen numbers due to its limitation of not being able to predict Knudsen layers.

\section{\label{sec:location}Location of singularities}
Finding an optimal location of singularity points has been a widely-discussed issue in the implementation of the MFS \cite{FR2002, chen2008, Alves2009, CKL2016, CNYC2023}.
It is well established that a high accuracy in the MFS is accompanied  with a high condition number due to the collocation matrix being ill-conditioned \citep{CH2020,
DML2009,CNYC2023,
WL2011}.
This is justified as the traditional condition number is not adequate for measuring the stability of the resulting system since the condition number does not take boundary data into account.
For instance, while forming matrix system \eqref{collocation_system}, the boundary data, such as $p_s$, $\tau_s$ or $\Delta \tau$, for the problems investigated in this paper appear in the vector $\mathcal{B}$ but not in the collocation matrix $\mathcal{M}$.
Hence, the (usual) condition number of the matrix $\mathcal{M}$ is not an adequate parameter to gauge the sensitivity of the MFS toward the location of the source points.

A more accurate estimation of the sensitivity of the MFS toward the location of the source points can be made by the \emph{effective condition number}, which also takes the boundary data into account (through the right-hand side vector).
The concept of the effective condition number has been used by many authors to determine an optimal location of the singularity points by conjecturing a reciprocal relationship between the inaccuracy of the MFS and the effective condition number \cite{DML2009, WL2011, CNYC2023}.
For both the problems discussed in the above sections, we have used the same strategy to place the source points.
Further details on this strategy are as follows.

Using the singular value decomposition, $\mathcal{M}$ (having dimensions $n\times m$) can be decomposed as
$\mathcal{M}=U D V^{\mathsf{T}}$, where $U$ and $V$ are $n\times n$ and $m \times m$ orthogonal matrices and $D$ is a $n\times m$ diagonal matrix containing the positive singular values in descending order: $\sigma_1 \geq \sigma_2 \geq \sigma_3 \geq \dots \geq \sigma_r>0$, where $r \leq m$.
The definitions of the (traditional) condition number $\kappa$ and the effective condition number $\kappa_{\mathrm{eff}}$ in $\ell^2$-norm are given by
\begin{align*}
\kappa=\frac{\sigma_1}{\sigma_r}\quad\text{and} \quad \kappa_{\mathrm{eff}}=\frac{\Vert\bm{\mathcal{B}} \Vert}{\sigma_r \Vert\bm{\mathcal{X}}\Vert}.
\end{align*}
Using the definition of the effective condition number,
we first verify the inverse relationship between the maximum error and the effective condition number.
\begin{figure}[!t]
\centering
\includegraphics[width=\columnwidth]{figures/err_eff_cond}
\caption{\label{error_keff}
The maximum absolute error $\epsilon_{\max}$ in the temperature and effective condition number $\kappa_{\mathrm{eff}}$ for the problem of flow between coaxial cylinders plotted over the dilation parameter $\alpha$ for $\mathrm{Kn}\approx 0$ and $N_b=N_s=100$.}
\end{figure}
Let $\alpha>1$ be the dilation parameter that determines the separation between the actual boundary (containing boundary nodes) and the fictitious boundary (containing singularities) such that $s_1=r_1/\alpha$ and $s_2=\alpha \,r_2$.
A larger value of $\alpha$ corresponds to a larger gap between the boundary nodes and source points.

For the problem described in Sec.~\ref{sec:prob1},
the maximum absolute error $\epsilon_{\max}$ in the temperature computed with the MFS and with the analytic solution for $\mathrm{Kn}\approx 0$ along with the effective condition number is plotted against the dilation parameter $\alpha$ in Fig.~\ref{error_keff}.
The figure shows that the inaccuracy of the MFS is roughly inversely proportional to the effective condition number.
It is also evident from the figure that the maximum value of the effective condition number is attained for $\alpha$ around $1.6$, where the effective condition number is of order $10^{8}$ and the absolute error is minimum.
It is worthwhile noting that the order of the effective condition number remains $10^{8}$ for higher values of $\alpha$ beyond $\alpha \approx 1.6$; similarly, the order of the maximum absolute error remains $10^{-5}$ for higher values of $\alpha$ beyond $\alpha \approx 1.6$.

In order to have further insight, the effective condition number for the problems considered in Sec.~\ref{sec:prob1} and Sec.~\ref{sec:prob2} is plotted against the dilation parameter in Fig.~\ref{keff} for different values of the Knudsen number.
The number of boundary nodes at either of the actual boundaries $N_b$ and the number of singularity points at either of the fictitious boundaries $N_s$ are taken as 100 (i.e.~$N_b=N_s=100$) in Fig.~\ref{keff}.
\begin{figure}[!tb]
\centering
\includegraphics[width=\columnwidth]{figures/eff_cond_coaxial.eps}
\\[3ex]
\includegraphics[width=\columnwidth]{figures/eff_cond_noncoaxial.eps}
\caption{\label{keff}
Variation of the effective condition number $\kappa_{\mathrm{eff}}$ with respect to the dilation parameter $\alpha$ for the problems described in Sec.~\ref{sec:prob1} (coaxial case) and Sec.~\ref{sec:prob2} (non-coaxial case).
The number of boundary nodes at either of the actual boundaries and the number of singularity points at either of the fictitious boundaries  are 100 (i.e.~$N_b=N_s=100$).}
\end{figure}
It can be noticed from the figure that the highest value of the effective condition number for a given Knudsen number is attained at a value of $\alpha$ somewhere in between $1.8$ and $2$.
Although, not shown here for succinctness, it turns out that the value of $\alpha$ at which the highest effective condition number is attained increases (decreases) with decrease (increase) in the number of boundary nodes and singularities.
Therefore, to save computational time, one can use smaller number of boundary nodes and source points along with a bigger value of $\alpha$.
From Fig.~\ref{keff}, although the effective condition number decreases on increasing $\alpha$ after a certain value of $\alpha$, we have not encountered any significant change in the results
on keeping the singularities farther  (or on taking bigger values of $\alpha$).
Therefore, it is apparently sufficient to just ensure $\alpha \geq 2$ to attain an optimal accuracy in the case of $N_b=N_s=100$.
Therefore, the fictitious boundaries for both problems have safely been positioned at locations for which $\alpha = 2$.

\section{\label{sec:conclusion}Conclusion and outlook}
The fundamental solutions of the CCR model in 2D have been determined by exploiting the fundamental solutions of some well-known partial differential equations, e.g.~the Laplace and biharmonic equations.
It turns out that the fundamental solutions of the linearized NSF and G13 equations in 2D can be recovered from the fundamental solutions of the CCR model in 2D derived in this paper by taking the coupling coefficient $\alpha_0$ as $0$ and $2/5$, respectively, in them.
The derived fundamental solutions for the two-dimensional CCR model have then been implemented in a numerical framework.

To gauge the capability of the developed numerical framework, two problems: (i) evaporating/condensing vapor flow between two coaxial cylinders, and (ii) temperature-driven rarefied gas flow between two non-coaxial cylinders having different temperatures, have been revisited.
These problems have already been investigated with the linearized BGK model in Refs.~[\onlinecite{Onishi1977} and \onlinecite{Aoki_Sone_1989}].
Comparing the results obtained from the MFS for the first problem with those from Refs.~[\onlinecite{Onishi1977}], the accuracy of the MFS with the CCR model in investigating rarefied gas flows with phase change is vivid, particularly for small Knudsen numbers.
Similarly, for the second problem, the results for the local flow fields, such as the temperature and velocity, obtained using the MFS with the CCR model compares quite well with those obtained using the linearized BGK model in the case of small Knudsen number; but the results for the local flow fields from the two models differ noticeably for larger Knudsen numbers, although their trends from both methods are qualitatively similar.
On the other hand, the MFS with the CCR model is able to capture the global flow fields, such as the drag force, quite accurately even for large Knudsen numbers.
In addition, since the MFS does not involve numerical computation of integrals and its implementation does not require the discretization of the domain, it is computationally efficient in comparison to the other numerical methods used for investigating rarefied gas flows.
This makes the MFS with the CCR model a favorable choice for investigating rarefied gas flows.
It should, however, be noted that  the position of singularity points plays a major role to achieve the best results.
By performing, effective condition number based studies for both problems, it has been established that the singularity points should be kept sufficiently far from the boundary nodes.

The utility of the derived fundamental solutions (and their implementation) can be perceived particularly for problems wherein the two-dimensional version of a problem is sufficient to study the complete problem in 3D (due to some symmetry) as the fundamental solutions of a model in 2D cannot be deduced directly from its counterpart in 3D and vice-versa.
Furthermore, the derived solutions can also be extended from monatomic to polyatomic gases (of Maxwell molecules) by taking $c_p=(5+\mathsf{n})/2$ and $\alpha_0=2/(5+\mathsf{n})$, and by choosing an appropriate value of the parameter $\mathsf{n}$ that denotes the degree of freedom in a polyatomic gas.

 Although the MFS is equally efficient for external flow problems as well (as demonstrated in Refs.~[\onlinecite{LC2016} and \onlinecite{RSCLS2021}]), we have not considered external flow problems in 2D as they are somewhat more involved in comparison to their 3D counterparts.
This is due to Stokes' paradox in which the solution diverges because of the presence of logarithmic term(s) in the two-dimensional fundamental solutions of Stokes' equation.
A similar logarithmic term appears in the two-dimensional fundamental solutions of the CCR model as well that makes the study of two-dimensional external flows with the MFS applied on the CCR model involved.
At present, we do not have a clear understanding of dealing with Stokes' paradox using the MFS for external flow problems, specifically for flows past an arbitrary geometry.
Notwithstanding, external flows with the MFS applied on the CCR model in 2D will be considered elsewhere in the future.

\section{Acknowledgement}
Himanshi gratefully acknowledges the financial support from the Council of Scientific and Industrial Research (CSIR) [File No.: 09/1022(0111)/2020-EMR-I]. A.S.R. acknowledges the financial support from the Science and Engineering Research Board, India through the grants SRG/2021/000790 and MTR/2021/000417.
Himanshi and V.K.G. also acknowledge Bhaskaracharya Mathematics Laboratory and Brahmagupta Mathematics Library supported by DST-FIST Project SR/FST/MS I/2018/26.