Current address: ]Department of Materials Science and Engineering, Cornell
University, Ithaca, New York 14853, USA

# Unraveling nano-scale effects of topotactic reduction in LaNiO2 crystals

Yu-Mi Wu<EMAIL_ADDRESS>[  Pascal Puphal  Masahiko Isobe  Bernhard
Keimer  Matthias Hepting  Y. Eren Suyolcu<EMAIL_ADDRESS>Peter A.
van Aken<EMAIL_ADDRESS>and<EMAIL_ADDRESS>Max Planck
Institute for Solid State Research, Heisenbergstrasse 1, 70569 Stuttgart,
Germany

###### Abstract

Infinite-layer nickelates stand as a promising frontier in the exploration of
unconventional superconductivity. Their synthesis through topotactic oxygen
reduction from the parent perovskite phase remains a complex and elusive
process. This study delves into the nano-scale effects of the topotactic
lattice transformation within LaNiO2 crystals. Leveraging high-resolution
scanning transmission electron microscopy and spectroscopy, our investigations
uncover a panorama of structural alterations, including grain boundaries and
coherent twin boundaries, triggered by reduction-induced transformations. In
addition, our analyses unveil the formation of an oxygen-rich disordered
transition phase encircling impurities and pervading crystalline domains and
the internal strain is accommodated by grain boundary formation. By unraveling
these nano-scale effects, our findings provide insights into the microscopic
intricacies of the topotactic reduction process elucidating the transition
from the perovskite to the infinite-layer phase within nickelate bulk
crystals.

††preprint: AIP/123-QED

##  I Introduction

The discovery of superconductivity in thin films of the infinite-layer (IL)
nickelate Nd0.8Sr0.2NiO2 Li  _et al._ (2019) has opened a new avenue to
explore a material system closely related to cuprate high-temperature
superconductors Keimer  _et al._ (2015). The structure of IL nickelates
contains quasi-two-dimensional NiO2 planes, obtainable through topotactic
oxygen reduction of the perovskite phase using either H2 gas or metal hydrides
as reducing agents Li  _et al._ (2019); Zeng  _et al._ (2020); Gao  _et al._
(2021); Hayward  _et al._ (1999); Crespin  _et al._ (2005). As an alternative
route for oxygen deintercalation, the deposition of an aluminum capping layer
on nickelate thin films has been employed recently Wei  _et al._ (2023a, b).
The nominal electronic configuration ($3d^{9}$) of the monovalent Ni1+ cations
in the NiO2 planes is isolelectronic to that of Cu2+ ions in undoped cuprates.
However, the electronic structure of the material appears to be different from
that of cuprates Hepting, Dean, and Lee (2021); Goodge  _et al._ (2021); Rossi
_et al._ (2021); Botana and Norman (2020); Been  _et al._ (2021); Held  _et
al._ (2022); Chen  _et al._ (2022).

The manifestation of superconductivity in IL nickelate thin films has
exhibited consistent reproducibility, extending recently across an entire
family of nickelates encompassing different rare-earth ions, $R$ = La and Pr
Osada  _et al._ (2020, 2021); Wang  _et al._ (2023), diverse levels of Sr and
Ca-substitution levels Li  _et al._ (2020a); Zeng  _et al._ (2022), and a
range of substrates Ren  _et al._ (2023); Lee  _et al._ (2023). Intriguingly,
superconductivity has also been observed in the quintuple-layer nickelate
Nd6Ni5O12 without Sr or Ca substitution, which can be considered as a self-
doped version of IL nickelates Pan  _et al._ (2022). Additionally, La3Ni2O7
was found to become superconducting under high pressure Sun  _et al._ (2023a),
although its $3d^{7.5}$ electronic configuration deviates significantly from
that of cuprates. This intricate panorama unveils tantalizing commonalities
and discrepancies between nickelates and cuprates, underscoring the need for a
comprehensive exploration to elucidate the full landscape of their behaviors.

In experimental practice, while the synthesis and fabrication of IL nickelates
present significant challenges, these processes have been extensively improved
for thin film samples Lee  _et al._ (2023, 2020); Sun  _et al._ (2023b); Chow
and Ariando (2022). Notably, during the topotactic reduction of thin films,
the underlying substrate provides epitaxial support essential for fostering
the IL phase formation. The thin film geometry, crucially, tends to yield a
single orientation of the reduced phase and minimizes the emergence of
impurity phases. Nevertheless, a considerable expansion of the in-plane
lattice constants and contraction of the out-of-plane lattice during the
reduction poses a challenge for epitaxial films. This phenomenon often prompts
the occurrence of extended defects as a means of alleviating epitaxial strain.
In particular, strain-relieving Ruddlesden-Popper type stacking faults often
occur densely in films grown on widely used SrTiO3 substrates Li  _et al._
(2020a); Osada  _et al._ (2021); Zeng  _et al._ (2022); Lee  _et al._ (2020);
Parzyck  _et al._ (2024). Other phases such as LaNiO2.5 and $a$-axis oriented
LaNiO2 can also form in response to reduction-induced strain Kawai  _et al._
(2010); Osada  _et al._ (2021). In general, the crystalline quality of the IL
structure is found to be closely linked to that of the perovskite precursor.
Furthermore, for uncapped films, a decomposition of the infinite-layer phase
in the uppermost regions of the films has been observed Lee  _et al._ (2020);
Onozuka  _et al._ (2016); Ikeda, Manabe, and Naito (2013), preventing surface-
sensitive techniques from the exploration of the structure-property
relationships within these materials.

As part of alternative approaches, recent efforts have directed attention
toward the synthesis of polycrystalline powders Li  _et al._ (2020b); Wang
_et al._ (2020); Ortiz  _et al._ (2022); Puphal  _et al._ (2022) and IL
nickelate single crystals Puphal  _et al._ (2021, 2023a). These single
crystals could, in principle, provide large single domains and cleavable
surfaces. However, in contrast to thin films benefiting from epitaxial
support, the arrangement of NiO2 planes within a crystal might align along any
of three equivalent pseudocubic [001] directions of the perovskite phase.
Consistent with this notion, three twin domains of the tetragonal $P4/mmm$
crystal structure were observed in topotactic La1-xCaxNiO2 crystals Puphal
_et al._ (2021). Recently, the synthesis of topotactic LaNiO2 crystals
exceeding one millimeter in size was achieved, utilizing a direct contact
method with the reducing agent CaH2 Puphal  _et al._ (2023a). This method was
subsequently optimized, employing cube-shaped LaNiO3 single crystals, each
with 1 mm edges, facilitating reproducible adjustment of the reduction
parameters Hayashida  _et al._ (2024). Despite electron backscatter
diffraction (EBSD) characterization of the resulting LaNiO2 crystals has
elucidated the distribution of 1-50 $\mu$m sized $P4/mmm$ twins across the
crystal surface Hayashida  _et al._ (2024), crucial insights into the
topotactically transformed microstructure in the bulk of these crystals remain
elusive. For a comprehensive understanding of the physical properties of
topotactically transformed systems, it is critical to gain insights into the
micro- and nanoscale impacts of the transformation process on the crystal
lattice, as the invasive nature of the process can result not only in
incompletely reduced regions within a crystal Puphal  _et al._ (2021), but
also in the emergence of unforeseen defects or phase decompositions Hu  _et
al._ (2024).

In the quest to unravel the complexities of the reduction process, spatially
resolved scanning tunneling electron microscopy (STEM) combined with electron
energy-loss spectroscopy (EELS) stands as a powerful tool to probe the local
lattice structure, elemental composition, and electronic structure down to the
atomic level, allowing for the identification and disentanglement of various
phases in topotactic systems Wu  _et al._ (2023). Here, we utilize the
capabilities of high-resolution STEM-EELS to scrutinize the structural effects
of the topotactic reduction in LaNiO2 crystals. Nanoscopic and microscopic
signatures of structural disorder and impurity phases throughout the reduced
crystal are unveiled by atomic-resolution STEM imaging. Our observations
demonstrate coherent twin boundaries with a relative rotation of $\sim$80∘ as
a consequence of strain accommodation from the pseudocubic to tetragonal phase
transformation. In addition, grain boundaries (GBs) formed by ultra-thin
La2O3-δ slabs and Ni layers in parallel to the main crystallographic
directions emerge within the crystal, likely originating from excessive
reduction. At a macroscopic scale, an intricate network of cracks and GB-like
regions separate single crystalline domains, while a structurally disordered
shell encircles the impurity phase, i.e. Ni inclusions. Our local
spectroscopic analyses reveal the oxygen-rich nature of these regions,
highlighting the formation of a LaNiO2+δ transition phase. Together, our
integrated structural and elemental analyses help to unveil the origin of
these defects within topotactically reduced crystals, illuminating the
intricate interplay between structure and strain in these systems.

##  II Methods

Single crystals of perovskite LaNiO3 were synthesized in a high-pressure
optical floating zone (OFZ) furnace (model HKZ, SciDre GmbH, Dresden,
Germany), as described in Refs. Puphal  _et al._ (2023a, b). The obtained
centimeter-sized crystals were cut into smaller cube-like shaped crystals,
with each surface of a cube corresponding to a pseudocubic [001] plane of the
rhombohedral $R\bar{3}c$ crystal structure of LaNiO3. The resulting cubes with
dimensions of $\sim$1 mm3 were subjected to a direct contact topotactic
reduction with CaH2, as described in Ref. Hayashida  _et al._ (2024).

Electron-transparent TEM specimens of the sample were prepared on a Thermo
Fisher Scientific focused ion beam (FIB) using the standard liftout method.
Samples with a size of 20 $\times$ 5 $\mu$m2 were thinned to 30 nm with 2 kV
Ga ions, followed by a final polish at 1 kV to reduce the effects of surface
damage. HAADF, ABF and EELS were recorded by a probe aberration-corrected JEOL
JEM-ARM200F scanning transmission electron microscope equipped with a cold-
field emission electron source and a probe Cs corrector (DCOR, CEOS GmbH), and
a Gatan K2 direct electron detector was used at 200 kV. STEM imaging and EELS
analyses were performed at probe semiconvergence angles of 20 and 28 mrad,
resulting in probe sizes of 0.8 and 1.0 Å, respectively. Collection angles for
STEM-HAADF and ABF images were 75 to 310 and 11 to 23 mrad, respectively. To
improve the signal-to-noise ratio of the STEM-HAADF and ABF data while
minimizing sample damage, a high-speed time series was recorded (2 $\mu$s per
pixel) and was then aligned and summed.

##  III Results

Figure 1: (a,c) SEM-secondary electron (SE) image of the perovskite LaNiO3 and
reduced LaNiO2 crystals, respectively. (b,d) Enlarged views of the images from
the orange and blue squares in (a) and (c), respectively. Inset zooms in the
region from the yellow square in (b), showing NiO impurities in the crystal.
(e) Atomic-resolution STEM-HAADF image taken in the LaNiO2 crystal after
topotactic reduction. The inset shows the corresponding fast Fourier transform
pattern. (f),(g) The simultaneously acquired STEM-HAADF and ABF images of the
region from the red rectangle in (e). (h)-(j) STEM-EELS elemental maps of La,
Ni, and the overlaid map, obtained from the green rectangle in (e).

Figures 1(a)-1(d) present scanning electron microscopy (SEM)-backscattered
electron (BSE) images of the crystals before and after the topotactic
reduction. A top-down view of a LaNiO3 crystal with typical lateral dimensions
of around 1 mm shows a polished surface without any visible domain boundaries
[Fig. 1(a)]. Upon closer inspection within a zoomed-in area, the emergence of
an irregularly shaped impurity phase, characterized by a darker contrast owing
to the presence of NiO impurities within the crystal, becomes evident [Fig.
1(b)] Puphal  _et al._ (2023b). While a LaNiO2 crystal with a minimal amount
of NiO inclusions was selected for the study in Ref. Hayashida  _et al._
(2024), for the present study, we chose a crystal with a high impurity
density. Figures 1(c) and 1(d) show the representative SEM-BSE images of a
LaNiO2 single crystal. Discernible impurity phases are observed within various
regions of this reduced crystal [Fig. 1(d)]. Furthermore, several adjacent
cracks and dark lines oriented in a parallel or orthogonal direction are
present, with a distance and a length ranging from a few to ten micrometers,
similar to the pattern previously observed in La1-xCaxNiO2 crystals Puphal
_et al._ (2021). These cracks and lines on the crystal surface correspond to
the dark amorphous GB-like regions with a width of a few hundred nanometers,
partitioning LaNiO2 domains within a cross-sectional TEM specimen (see Fig. S1
of the Supplemental Material for more details). The existence of these cracks
after reduction indicates a consequence of topotactically induced structural
changes.

To understand the local structural changes upon reduction in more detail, we
first investigate the atomic lattice in the LaNiO2 crystal by atomic-
resolution STEM high-angle annular dark-field (HAADF) imaging. Identical TEM
specimens from both impurity-free and impurity-containing regions of LaNiO2
were prepared for comparison. Within a single-crystalline domain in the region
without any impurities, the STEM-HAADF image presents a high crystalline
quality without any defects or impurity phase after the reduction process
[Fig. 1(e)]. Furthermore, images with a larger field of view do not reveal any
regions with defects or impurity phases in one individual domain (Fig. S1).
The typical domain size within the crystal is found to be at a scale of a few
tenths of nanometers. The crystal symmetry is shown by the corresponding fast
Fourier transforms (FFT) of the image [Fig. 1(e) inset]. The difference of the
distances between the FFT maxima in reciprocal space along the [001] and [100]
axis elucidates a tetragonality $c/a$ ratio of approximately 0.8 (highlighted
by yellow arrows on FFT patterns). This indicates a tetragonal structure with
an out-of-plane lattice contraction and an in-plane lattice expansion upon
topotactic removal of apical oxygen atoms. This $c/a$ ratio aligns with the
lattice parameters $a,b=3.9642(3)$ Å and $c=3.3561(3)$ Å determined by x-ray
diffraction for a LaNiO2 crystal in Ref. Hayashida  _et al._ (2024). Further
zooming-in STEM-HAADF and annular bright-field (ABF) images are displayed in
Fig. 1(f) and 1(g). The distribution of oxygen ions including filled and empty
apical oxygen sites is clearly visible by ABF imaging. The absence of image
contrast at apical oxygen sites confirms the infinite-layer structure of
LaNiO2 crystals. EELS elemental maps of La and Ni obtained from the crystal
are shown in Figs. 1(h)-1(j) using La $M_{5,4}$ and Ni $L_{3,2}$ edges,
respectively. The maps confirm the uniform distributions of La and Ni
throughout the crystal structure. This suggests that each domain retains a
stoichiometric and high-quality infinite-layer phase within each domain
following the reduction process.

Figure 2: (a) A low-magnification STEM-HAADF image of a LaNiO2 crystal. The
blue and red dashed lines define two types of boundaries observed in the
crystal. (b) The zoom-in image of the area from the blue square in (a). The
insets are the fast Fourier transform patterns of grains 1 and 2, showing the
same structural phases between the boundary. (c) Enlarged view of the area
from the yellow rectangle in (b) and a sketch of the structural model for the
coherent twin boundary. (d) Enlarged view of the area from the red rectangle
in (a). The green and red shaded areas indicate different atomic stacking from
the single-crystalline regions. (e)-(g) STEM-EELS elemental maps of La, Ni, O
from the yellow rectangle in (d) across the grain boundary. (h) Sketch of the
structural model of the grain boundary formed by Ni metal and La2O3-δ layers.

The alternation in the tetragonality $c/a$ ratio of the crystal through
reduction can destabilize the infinite-layer structure in LaNiO2, potentially
causing structural inhomogeneities or defects Puphal  _et al._ (2023a). Hence,
our focus shifts toward comprehending the structural changes in the crystal
after the reduction process. At a macroscopic scale, the crystal presents a
network of GB-like regions arranged orthogonally, suggesting that LaNiO2
retains aligned domains with $a$ or $c$-axis orientations globally through
reduction (Fig. S1). This implies an out-of-plane (in-plane) lattice
contraction (expansion) of the crystal upon reduction paralleling any of the
three symmetry-related [100] axes of the pseudocubic perovskite phase. Zooming
into a low-magnification STEM-HAADF image of LaNiO2 in Fig. 2(a), we observe
inclined and intersecting dark features. The dashed lines denote two types of
boundaries: lines oriented diagonally (blue); and lines oriented vertically
and horizontally (red). Magnification of an area along one of the dark lines
[Fig. 2(b)], as denoted by the blue dashed line and arrows, reveals two grains
oriented along [100] and [001] directions [Fig. 2(b), insets]. These two
crystalline grains, labeled 1 and 2, exhibit an identical lattice structure
but with a rotation relative to one another around the misorientation axis. In
a detailed view of the boundary between two grains in Figure 2(c), grain 2
appears rotated $\sim$80∘ around the [001] axis with respect to grain 1. Such
twinning, where the crystal structure of one part mirrors the matrix, often
emerges during a phase transformation within a material to accommodate the
strain Christian and Mahajan (1995). In LaNiO2 crystals, this arises from the
perovskite pseudocubic to a tetragonal phase transition induced by topotactic
reduction. Using simple geometrical calculations, the rotation angle,
$\theta$, is determined as: $\theta=2\tan^{-1}(c/a)$. The observed $c/a$ ratio
of $\sim$0.8 from the FFT maxima in Fig. 1(e) aligns with a rotation angle of
80∘ between the matrix around the [001] axis. This is consistent with the
twinning angle observed from the image and confirms that the pseudocubic to
tetragonal phase transformation resulting in an internal strain, causes
twinning throughout the crystal. The dark contrast of the boundary can be
attributed to diffuse and dynamic scattering resulting from the de-channeling
of the electron beam as a consequence of the difference in the structure of
the boundary compared to the bulk domains Wang (1994); Yu, Muller, and Silcox
(2004).

Figure 3: (a) Low-magnification STEM-HAADF image showing the NiO impurity in
the LaNiO3 crystal. The inset shows a magnified image of the interface region
between LaNiO3 and the impurity from the green rectangle of the blue rectangle
in (a). (b) Low-magnification STEM-HAADF image showing the Ni inclusion
surrounded by a shell in the LaNiO2 crystal. The right panel shows the
magnified view of the Ni inclusion. (c)-(f) STEM-EELS elemental maps of La,
Ni, O and the superimposed map obtained from the orange rectangle in (b),
including crystalline and grain boundary-like regions in LaNiO2, oxygen-rich
shell and Ni.

On the other hand, as indicated by the red dashed lines in Fig. 2(a), a
different type of boundary becomes evident, aligned parallel to
crystallographic directions. As shown in Fig. 2(d), the dark line aligns along
the [100] direction of LaNiO2. Contrary to the grains connected through a twin
boundary as shown in Fig. 2(b)-2(c), grains here share an interfacial
transition region structurally and chemically distinct from their adjacent
counterparts. Specifically, a Ni layer with a cubic rock-salt structure is
sandwiched between additional La-O layers exhibiting a horizontal relative
shift by $a$/2[100]. At the GB region, the spacing between the La-O layer and
its adjacent one ($\sim$2.8 Å) is smaller than the spacing between two
neighboring La-O layers ($\sim$3.4 Å) in the bulk LaNiO2 phase. To confirm the
chemical composition, STEM-EELS elemental maps of La, Ni, and O across the GB
are displayed in Fig. 2(e)-2(g). The La map reveals the extra LaO1-δ layers
segregating the LaNiO2 from the inner Ni layer. The region with enhanced Ni
content indicates the presence of the Ni layer at the GB, while the oxygen
intensity decreases at the GB region. Previous powder x-ray diffraction
measurements Ref. Puphal  _et al._ (2023a) have identified a chemical
decomposition of LaNiO2 crystals into elemental Ni, La2O3, and LaOH due to
excessive reduction. Consequently, our findings likely reveal the
decomposition process associated with the formation of GB containing elemental
Ni bordered by La2O3-δ blocks between domains.

Having discussed the local structural deformations, we now turn to the larger-
scale imperfections in the crystal due to topotactic reduction. As shown in
Fig. 1(b) and 1(d), impurity phases are present in both unreduced and reduced
crystals. To discern the origin of these impurities, we first investigate the
unreduced crystal. The low-magnification STEM-HAADF image of LaNiO3 in Fig.
3(a) reveals a pristine crystalline region housing a rod-shaped NiO impurity.
Upon closer inspection, the atomic structure of cubic rock-salt NiO is in
close contact with pseudocubic perovskite LaNiO3 with a structurally coherent
interface due to closely aligned lattice symmetries [Fig. 3(a), inset].
Supporting STEM-HAADF and EELS maps in Fig. S2 confirm the NiO impurity phase
and a chemically sharp LaNiO3/NiO interface. After reduction, an extent of
decomposition near the impurity and within the crystalline region throughout
the crystal appears, while a partially polycrystalline character is present
inside the Ni inclusion [Fig. 3(b)]. Zooming in on a single grain exhibits a
cubic rock-salt structure within the impurity phase. Interestingly, a dark
shell with a width of approximately a hundred nanometer surrounding the Ni
inclusion separates the impurity phase from the crystalline region, unlike the
sharp LaNiO3/NiO interface observed pre-reduction. Within the LaNiO2 bulk
region, some vertical lines and irregular dark GB-like regions with a width of
a few ten nanometers are observed. Both the shell and GB-like region exhibit
an amorphous structure, displaying dark contrast in the image owing to the
diffuse scattering Wang (1994).

In Fig. 3(c)-3(f), the elemental maps reveal how the crystal structure is
modified due to the reduction process: a Ni-rich phase emerges without other
contributions confirming that the NiO impurity in LaNiO3 is reduced to Ni
metal in LaNiO2. An excess of Ni accompanying La and O deficiencies are also
present at the vertical lines in the bulk region, indicating a formation of Ni
metal layers that are identical to the structural modifications observed upon
excess reduction (c.f., Fig. 2(d)). Besides, the O elemental map in Fig. 3(e)
provides a clear view of enhanced intensity in the amorphous shell and GB-like
regions. These findings suggest that the removal of oxygen drives the
transformation of NiO impurities into polycrystalline Ni inclusions. The
observed structural disorder within the surrounding shell and GB-like region
blocks oxygen to be pulled out in the crystal during the reduction and,
therefore appears to be oxygen-rich relative to stoichiometric LaNiO2 domains,
leading to the formation of an amorphous LaNiO2+δ phase (Fig. S3). Such a
disordered phase can potentially be associated with a mechanism relieving
reduction-induced strain as well. Functioning as a transition layer, the
LaNiO2+δ phase between the domains and Ni impurities and in the crystal relax
the local strain arising from the lattice mismatch between the defect and bulk
regions.

##  IV Conclusion

In summary, we scrutinize how the topotactic oxygen reduction impacts the
structural transformation of LaNiO3 to the infinite-layer phase. Formation of
macroscopic and local structural imperfections are observed in response to the
reduction-induced crystal symmetry changes in LaNiO2. At the nanometer scale,
we discern two types of boundaries separating domains. Globally, the infinite-
layer phase is preserved in LaNiO2, featuring coherent twin boundaries
oriented diagonally to accommodate the internal strain. Besides, non-
stoichiometric grain boundaries coexist within the crystal, alongside a
decomposition into La2O3-δ and Ni metal due to excessive reduction, formed in
parallel to crystallographic orientations between crystalline regions. Moving
to the micrometer scale, reduction-induced structural disorder surrounds the
Ni metal impurities and resides within crystalline domains. These disordered
regions, impeding the local oxygen reduction, result in a transition LaNiO2+δ
phase. Our atomic-scale observation of these defects provides a good basis for
comprehending the reduction mechanism and offers valuable insights into
optimizing the synthesis of nickelate single crystals and related materials in
future endeavors.

###### Acknowledgements.

We thank J. Deuschle for TEM specimen preparation. This project has received
funding in part from the European Union’s Horizon 2020 research and innovation
program under Grant Agreement No. 823717–ESTEEM3.

## References

  * Li  _et al._ (2019) D. Li, K. Lee, B. Y. Wang, M. Osada, S. Crossley, H. R. Lee, Y. Cui, Y. Hikita, and H. Y. Hwang, “Superconductivity in an infinite-layer nickelate,” Nature 572, 624–627 (2019).
  * Keimer  _et al._ (2015) B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, “From quantum matter to high-temperature superconductivity in copper oxides,” Nature 518, 179–186 (2015).
  * Zeng  _et al._ (2020) S. Zeng, C. S. Tang, X. Yin, C. Li, M. Li, Z. Huang, J. Hu, W. Liu, G. J. Omar, H. Jani, _et al._ , “Phase diagram and superconducting dome of infinite-layer Nd1-xSrxNiO2 thin films,” Phys. Rev. Lett. 125, 147003 (2020).
  * Gao  _et al._ (2021) Q. Gao, Y. Zhao, X.-J. Zhou, and Z. Zhu, “Preparation of Superconducting Thin Films of Infinite-Layer Nickelate Nd0.8Sr0.2NiO2,” Chinese Physics Letters 38, 077401 (2021).
  * Hayward  _et al._ (1999) M. A. Hayward, M. A. Green, M. J. Rosseinsky, and J. Sloan, “Sodium Hydride as a Powerful Reducing Agent for Topotactic Oxide Deintercalation Synthesis and Characterization of the Nickel(I) Oxide LaNiO2,” J. Am. Chem. Soc. 121, 8843–8854 (1999).
  * Crespin  _et al._ (2005) M. Crespin, O. Isnard, F. Dubois, J. Choisnet, and P. Odier, “LaNiO2: Synthesis and structural characterization,” J. Solid State Chem. 178, 1326 – 1334 (2005).
  * Wei  _et al._ (2023a) W. Wei, K. Shin, H. Hong, Y. Shin, A. S. Thind, Y. Yang, R. F. Klie, F. J. Walker, and C. H. Ahn, “Solid state reduction of nickelate thin films,” Phys. Rev. Mater. 7, 013802 (2023a).
  * Wei  _et al._ (2023b) W. Wei, D. Vu, Z. Zhang, F. J. Walker, and C. H. Ahn, “Superconducting Nd1-xEuxNiO2 thin films using in situ synthesis,” Sci. Adv. 9, eadh3327 (2023b).
  * Hepting, Dean, and Lee (2021) M. Hepting, M. P. M. Dean, and W.-S. Lee, “Soft X-Ray spectroscopy of low-valence nickelates,” Front. Phys. 9:808683 (2021).
  * Goodge  _et al._ (2021) B. H. Goodge, D. Li, K. Lee, M. Osada, B. Y. Wang, G. A. Sawatzky, H. Y. Hwang, and L. F. Kourkoutis, “Doping evolution of the Mott–Hubbard landscape in infinite-layer nickelates,” Proc. Natl. Acad. Sci. 118, e2007683118 (2021).
  * Rossi  _et al._ (2021) M. Rossi, H. Lu, A. Nag, D. Li, M. Osada, K. Lee, B. Y. Wang, S. Agrestini, M. Garcia-Fernandez, J. J. Kas, Y. D. Chuang, Z. X. Shen, H. Y. Hwang, B. Moritz, K. J. Zhou, T. P. Devereaux, and W. S. Lee, “Orbital and spin character of doped carriers in infinite-layer nickelates,” Phys. Rev. B 104, L220505 (2021).
  * Botana and Norman (2020) A. S. Botana and M. R. Norman, “Similarities and differences between LaNiO2 and CaCuO2 and implications for superconductivity,” Phys. Rev. X 10, 011024 (2020).
  * Been  _et al._ (2021) E. Been, W.-S. Lee, H. Y. Hwang, Y. Cui, J. Zaanen, T. Devereaux, B. Moritz, and C. Jia, “Electronic structure trends across the rare-earth series in superconducting infinite-layer nickelates,” Phys. Rev. X 11, 011050 (2021).
  * Held  _et al._ (2022) K. Held, L. Si, P. Worm, O. Janson, R. Arita, Z. Zhong, J. M. Tomczak, and M. Kitatani, “Phase diagram of nickelate superconductors calculated by dynamical vertex approximation,” Front. Phys. 9:810394 (2022).
  * Chen  _et al._ (2022) H. Chen, A. Hampel, J. Karp, F. Lechermann, and A. J. Millis, “Dynamical mean field studies of infinite layer nickelates: Physics results and methodological implications,” Front. Phys. 10:835942 (2022).
  * Osada  _et al._ (2020) M. Osada, B. Y. Wang, B. H. Goodge, K. Lee, H. Yoon, K. Sakuma, D. Li, M. Miura, L. F. Kourkoutis, and H. Y. Hwang, “A superconducting praseodymium nickelate with infinite layer structure,” Nano Lett. 20, 5735–5740 (2020).
  * Osada  _et al._ (2021) M. Osada, B. Y. Wang, B. H. Goodge, S. P. Harvey, K. Lee, D. Li, L. F. Kourkoutis, and H. Y. Hwang, “Nickelate Superconductivity without Rare-Earth Magnetism: (La,Sr)NiO2,” Adv. Mater. 33, 2104083 (2021).
  * Wang  _et al._ (2023) B. Y. Wang, T. C. Wang, Y.-T. Hsu, M. Osada, K. Lee, C. Jia, C. Duffy, D. Li, J. Fowlie, M. R. Beasley, T. P. Devereaux, I. R. Fisher, N. E. Hussey, and H. Y. Hwang, “Effects of rare-earth magnetism on the superconducting upper critical field in infinite-layer nickelates,” Sci. Adv 9, eadf6655 (2023).
  * Li  _et al._ (2020a) D. Li, B. Y. Wang, K. Lee, S. P. Harvey, M. Osada, B. H. Goodge, L. F. Kourkoutis, and H. Y. Hwang, “Superconducting dome in ${\mathrm{Nd}}_{1-x}{\mathrm{Sr}}_{x}{\mathrm{NiO}}_{2}$ infinite layer films,” Phys. Rev. Lett. 125, 027001 (2020a).
  * Zeng  _et al._ (2022) S. Zeng, C. Li, L. E. Chow, Y. Cao, Z. Zhang, C. S. Tang, X. Yin, Z. S. Lim, J. Hu, P. Yang, and A. Ariando, “Superconductivity in infinite-layer nickelate La1-xCaxNiO2 thin films,” Sci. Adv. 8, eabl9927 (2022).
  * Ren  _et al._ (2023) X. Ren, J. Li, W.-C. Chen, Q. Gao, J. J. Sanchez, J. Hales, H. Luo, F. Rodolakis, J. L. McChesney, T. Xiang, J. Hu, R. Comin, Y. Wang, X. Zhou, and Z. Zhu, “Possible strain-induced enhancement of the superconducting onset transition temperature in infinite-layer nickelates,” Commun. Phys. 6, 341 (2023).
  * Lee  _et al._ (2023) K. Lee, B. Y. Wang, M. Osada, B. H. Goodge, T. C. Wang, Y. Lee, S. Harvey, W. J. Kim, Y. Yu, C. Murthy, _et al._ , “Linear-in-temperature resistivity for optimally superconducting (Nd,Sr) NiO2,” Nature 619, 288–292 (2023).
  * Pan  _et al._ (2022) G. A. Pan, D. Ferenc Segedin, H. LaBollita, Q. Song, E. M. Nica, B. H. Goodge, A. T. Pierce, S. Doyle, S. Novakov, D. Córdova Carrizales, _et al._ , “Superconductivity in a quintuple-layer square-planar nickelate,” Nature Materials 21, 160–164 (2022).
  * Sun  _et al._ (2023a) H. Sun, M. Huo, X. Hu, J. Li, Z. Liu, Y. Han, L. Tang, Z. Mao, P. Yang, B. Wang, _et al._ , “Signatures of superconductivity near 80 k in a nickelate under high pressure,” Nature 621, 493–498 (2023a).
  * Lee  _et al._ (2020) K. Lee, B. H. Goodge, D. Li, M. Osada, B. Y. Wang, Y. Cui, L. F. Kourkoutis, and H. Y. Hwang, “Aspects of the synthesis of thin film superconducting infinite-layer nickelates,” Apl Materials 8 (2020).
  * Sun  _et al._ (2023b) W. Sun, Y. Li, R. Liu, J. Yang, J. Li, W. Wei, G. Jin, S. Yan, H. Sun, W. Guo, Z. Gu, Z. Zhu, Y. Sun, Z. Shi, Y. Deng, X. Wang, and Y. Nie, “Evidence for anisotropic superconductivity beyond Pauli limit in infinite-layer lanthanum nickelates,” Adv. Mater. 35, 2303400 (2023b).
  * Chow and Ariando (2022) L. E. Chow and A. Ariando, “Infinite-layer nickelate superconductors: A current experimental perspective of the crystal and electronic structures,” Front. Phys. 10:834658 (2022).
  * Parzyck  _et al._ (2024) C. T. Parzyck, V. Anil, Y. Wu, B. H. Goodge, M. Roddy, L. F. Kourkoutis, D. G. Schlom, and K. M. Shen, “Synthesis of thin film infinite-layer nickelates by atomic hydrogen reduction: Clarifying the role of the capping layer,” APL Materials 12, 031132 (2024).
  * Kawai  _et al._ (2010) M. Kawai, K. Matsumoto, N. Ichikawa, M. Mizumaki, O. Sakata, N. Kawamura, S. Kimura, and Y. Shimakawa, “Orientation change of an infinite-layer structure LaNiO2 epitaxial thin film by annealing with CaH2,” Crystal growth & design 10, 2044–2046 (2010).
  * Onozuka  _et al._ (2016) T. Onozuka, A. Chikamatsu, T. Katayama, T. Fukumura, and T. Hasegawa, “Formation of defect-fluorite structured NdNiOxHy epitaxial thin films via a soft chemical route from NdNiO3 precursors,” Dalton Transactions 45, 12114–12118 (2016).
  * Ikeda, Manabe, and Naito (2013) A. Ikeda, T. Manabe, and M. Naito, “Improved conductivity of infinite-layer LaNiO2 thin films by metal organic decomposition,” Physica C: Superconductivity 495, 134–140 (2013).
  * Li  _et al._ (2020b) Q. Li, C. He, J. Si, X. Zhu, Y. Zhang, and H.-H. Wen, “Absence of superconductivity in bulk Nd1-xSrxNiO2,” Commun. Mater. 1, 16 (2020b).
  * Wang  _et al._ (2020) B.-X. Wang, H. Zheng, E. Krivyakina, O. Chmaissem, P. P. Lopes, J. W. Lynn, L. C. Gallington, Y. Ren, S. Rosenkranz, J. F. Mitchell, and D. Phelan, “Synthesis and characterization of bulk ${\mathrm{Nd}}_{1-x}{\mathrm{Sr}}_{x}\mathrm{Ni}{\mathrm{O}}_{2}$ and ${\mathrm{Nd}}_{1-x}{\mathrm{Sr}}_{x}\mathrm{Ni}{\mathrm{O}}_{3}$,” Phys. Rev. Mater. 4, 084409 (2020).
  * Ortiz  _et al._ (2022) R. A. Ortiz, P. Puphal, M. Klett, F. Hotz, R. K. Kremer, H. Trepka, M. Hemmida, H.-A. K. von Nidda, M. Isobe, R. Khasanov, H. Luetkens, P. Hansmann, B. Keimer, T. Schäfer, and M. Hepting, “Magnetic correlations in infinite-layer nickelates: An experimental and theoretical multimethod study,” Phys. Rev. Research 4, 023093 (2022).
  * Puphal  _et al._ (2022) P. Puphal, V. Pomjakushin, R. A. Ortiz, S. Hammoud, M. Isobe, B. Keimer, and M. Hepting, “Investigation of hydrogen incorporations in bulk infinite-layer nickelates,” Front. Phys. 10:842578 (2022).
  * Puphal  _et al._ (2021) P. Puphal, Y.-M. Wu, K. Fuersich, H. Lee, M. Pakdaman, J. A. Bruin, J. Nuss, Y. E. Suyolcu, P. A. van Aken, B. Keimer, _et al._ , “Topotactic transformation of single crystals: From perovskite to infinite-layer nickelates,” Sci. Adv. 7, eabl8091 (2021).
  * Puphal  _et al._ (2023a) P. Puphal, B. Wehinger, J. Nuss, K. Küster, U. Starke, G. Garbarino, B. Keimer, M. Isobe, and M. Hepting, “Synthesis and physical properties of ${\mathrm{lanio}}_{2}$ crystals,” Phys. Rev. Mater. 7, 014804 (2023a).
  * Hayashida  _et al._ (2024) S. Hayashida, V. Sundaramurthy, P. Puphal, M. Garcia-Fernandez, K.-J. Zhou, B. Fenk, M. Isobe, M. Minola, Y.-M. Wu, Y. E. Suyolcu, P. A. van Aken, B. Keimer, and M. Hepting, “Investigation of spin excitations and charge order in bulk crystals of the infinite-layer nickelate ${\mathrm{LaNiO}}_{2}$,” Phys. Rev. B 109, 235106 (2024).
  * Hu  _et al._ (2024) K. Hu, Q. Li, D. Song, Y. Jia, Z. Liang, S. Wang, H. Du, H.-H. Wen, and B. Ge, “Atomic scale disorder and reconstruction in bulk infinite-layer nickelates lacking superconductivity,” Nat. Commun. 15, 5104 (2024).
  * Wu  _et al._ (2023) Y.-M. Wu, P. Puphal, H. Lee, J. Nuss, M. Isobe, B. Keimer, M. Hepting, Y. E. Suyolcu, and P. A. van Aken, “Topotactically induced oxygen vacancy order in nickelate single crystals,” Phys. Rev. Mater. 7, 053609 (2023).
  * Puphal  _et al._ (2023b) P. Puphal, V. Sundaramurthy, V. Zimmermann, K. Küster, U. Starke, M. Isobe, B. Keimer, and M. Hepting, “Phase formation in hole- and electron-doped rare-earth nickelate single crystals,” APL Materials 11, 081107 (2023b).
  * Christian and Mahajan (1995) J. W. Christian and S. Mahajan, “Deformation twinning,” Progress in Materials Science 39, 1–157 (1995).
  * Wang (1994) Z. Wang, “Dislocation contrast in high-angle hollow-cone dark-field TEM,” Ultramicroscopy 53, 73–90 (1994).
  * Yu, Muller, and Silcox (2004) Z. Yu, D. A. Muller, and J. Silcox, “Study of strain fields at a-Si/c-Si interface,” Journal of Applied Physics 95, 3362–3371 (2004).


# A velocity-space adaptive unified gas kinetic scheme for continuum and
rarefied flows

Tianbai Xiao<EMAIL_ADDRESS>Kun Xu<EMAIL_ADDRESS>Qingdong Cai
<EMAIL_ADDRESS>Department of Mechanics and Engineering Science, College of
Engineering, Peking University, Beijing 100871, China  Department of
Mathematics, Department of Mechanical and Aerospace Engineering, Hong Kong
University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong
Shenzhen Research Institute, Hong Kong University of Science and Technology,
Shenzhen 518057, China

###### Abstract

The compressible flow has intrinsically multiple scale nature due to the large
variations of gas density and characteristic scale of flow structure,
especially in hypersonic and reentry problems. It is challenging to construct
an accurate and efficient numerical algorithm to capture non-equilibrium flow
physics across different regimes. In this paper, a unified gas kinetic scheme
with adaptive velocity space (AUGKS) for multiscale flow transport will be
developed. In near-equilibrium flow regions, particle distribution function is
close to the Chapman-Enskog expansion and can be formulated with a continuous
velocity space, where only macroscopic conservative variables are updated.
With the emerging of non-equilibrium effects, the AUGKS automatically switches
to a discrete velocity space to follow the evolution of particle distribution
function. Based on the Chapman-Enskog expansion, a criterion is proposed in
this paper to quantify the intensity of non-equilibrium effects and is used
for the continuous-discrete velocity space transformation. Following the
scale-dependent local evolution solution, the AUGKS presents the discretized
gas dynamic equations directly on the cell size and time step scales, i.e.,
the so-called direct modeling method. As a result, the scheme is able to
capture the cross-scale flow physics from particle transport to hydrodynamic
wave propagation, and provides a continuous variation of solutions from the
Boltzmann to the Navier-Stokes. Under the unified framework, different from
conventional DSMC-NS hybrid method, the AUGKS does not need a buffer zone to
match up kinetic and hydrodynamic solutions. Instead, a continuous and
discrete particle velocity space is naturally connected, which is feasible for
the numerical simulations with unsteadiness or complex geometries. Compared
with the asymptotic preserving (AP) methods which solves kinetic equations
uniformly over the entire flow field with discretized velocity space, the
current velocity-space adaptive unified scheme speeds up the computation and
reduces the memory requirement in multiscale flow problems, and maintains the
equivalent accuracy. The AUGKS provides an effective tool for non-equilibrium
flow simulations.

###### keywords:

unified gas kinetic scheme, multiscale flow, non-equilibrium phenomena,
adaptive velocity space

##  1 Introduction

The gaseous flow shows a diverse set of behaviors on different characteristic
scales. For example, within the mean free path and collision time of gas
molecules, particles travel freely during most of time with rare
intermolecular collisions, leading to peculiar non-equilibrium flow dynamics.
However, on a much macroscopic level, the accumulating effect of collisions
results in an equalization of local temperature and velocity, where the
moderate non-equilibrium effects can be well described by viscous transport,
heat conduction and mass diffusion, i.e., the so called transport phenomena
[1]. From microscopic particle transport to macroscopic fluid motion, there is
a continuous variation of flow dynamics. Generally, different flow regimes can
be categorized qualitatively according to the Knudsen number, which is defined
as the ratio of the molecular mean free path to a characteristic physical
length scale. With the variation of $Kn$, the whole flow domain can be
qualitatively divided into continuum ($Kn<0.001$), slip ($0.001<Kn<0.1$),
transition ($0.1<Kn<10$) and free molecular regimes ($Kn>10$) [2]. When $Kn$
is large, the particle transport and collision can be depicted separately in
the Boltzmann equation. In another limit with extremely small $Kn$, the
Navier-Stokes equations are routinely used to describe macroscopic flow.

The traditional computational fluid dynamics targets to get numerical
solutions of the corresponding governing equations. For example, the most
widely used numerical methods for the Boltzmann equation are the direct
Boltzmann solvers [3] and the direct simulation Monte Carlo (DSMC) method [4].
In the former methodology, a discretized particle velocity space is
constructed and the particle distribution function is updated from transport
and collision terms respectively. On the other hand, the DSMC method mimics
the same physical process while the distribution function is now represented
by a large amount of test particles and the collision term is calculated
statistically. Due to the splitting treatment of particle transport and
collision, the mesh size and time step should be restricted by the mean free
path and collision time, and the computational cost is proportional to the
amount of discretized velocity points or test particles used in the
simulation. Meanwhile, the compressible Navier-Stokes solvers are mostly based
on the Riemann solvers for inviscid flux and the central difference method for
viscous terms. The macroscopic flow variables are followed in the simulation.
Compared with the kinetic methods, the computational cost of continuum flow
solvers is much reduced.

The rapid developments of aerospace industry face new challenges for accurate
and efficient simulation of complex flows. For example, when a shuttle
reenters into the atmosphere, the surrounding gas becomes denser and denser
from the rarefied upper atmosphere to lower continuum region, and thus the
vehicle sees cross-scale flow physics during the landing process. Besides,
localized non-equilibrium flow structures often emerge around the vehicle in
hypersonic cruise as a result of the geometric effect, such as shock,
rarefaction wave, boundary layer and wake turbulence. It is natural to couple
different numerical methods in different domains to calculate aerodynamic
force and heat efficiently. Therefore, hybrid algorithms which combine
continuum and kinetic approaches have been developed to simulate multiple
scale flows where continuum and rarefied flow physics coexist in a single flow
simulation [5, 6, 7, 8, 9, 10, 11, 12, 13]. In these numerical schemes, the
main flow structure of the flow field is simulated by the continuum methods
efficiently, with highly dissipative non-equilibrium regions being resolved by
the kinetic methods. Due to the complicated fivefold collision integral in the
Boltzmann equation, the prevailing kinetic solver used for hybridization
currently is the DSMC method. In the calculation, a dynamic parameter is
needed to determine the breakdown point of continuum description and decompose
the flow domain. It is easy to implement the hybrid method into a parallel
simulation since the physical domain has already been splitted into blocks on
different computational nodes.

In kinetic theory, the Chapman-Enskog expansion [14] bridges the NS and
Boltzmann solutions. Although this successive expansion is mathematically
attractive, there is little information provided about the intrinsic scale for
the validation of macroscopic equations, such as the specific fluid element in
the NS modeling. The success of the mathematical derivation of high-order
equations with inclusion of the so-called Burnett or super-Burnett terms is
limited due to the lack of specified modeling scales for these extended
hydrodynamic equations. Besides, due to the uncertainty in choosing the length
scale in the definition of Knudsen number, it becomes rather tough to predict
a universal breakdown criterion for the Chapman-Enskog expansion and the use
of the NS solutions, although it is defined empirically that the NS equations
are valid when ${Kn}\leq 0.01$. In addition, on particle mean free path and
collision time scales, the kinetic method provides much more degrees of
freedom to describe gas dynamics, while in the buffer zone these information
needs to be shrunk to the mean field variables on a macroscopic level, such as
density, momentum, energy, stress and heat flux. The inherent incompatibility
between the particle-based and PDE-based methods leads to considerable
difficulties in the hybridization. Usually a buffer zone is delicately defined
for the information exchange between kinetic and continuum solutions, which
may become more difficult in the simulation of unsteady flows.

In recent years, the unified gas kinetic scheme (UGKS) has been developed for
cross-scale flows [15, 16]. Based on the direct modeling on the mesh size
scale, the time-evolving flux function from the kinetic equation provides a
smooth transition from particle transport to hydrodynamic wave propagation.
The UGKS is an asymptotic preserving (AP) scheme, which preserves the discrete
analogy of the Chapman-Enskog expansion and the time step in the simulation is
not restricted by the collision time [17]. However, for the UGKS the memory
requirement and computational cost due to the discretized velocity space
limits its wide applications in aerospace industry. In this paper, we develop
an adaptive unified gas kinetic scheme (AUGKS), where the continuous and
discrete velocity space is transformed dynamically in a single framework. In
the near-equilibrium region, the Chapman-Enskog expansion is used for the
construction of distribution function with a continuous velocity space, and
thus only macroscopic conservative variables are stored and updated in the
simulation. With increasing non-equilibrium effects, the AUGKS tracks the
evolution of distribution function directly with a discrete velocity space.
Based on the Chapman-Enskog expansion, a criterion to switch continuous-
discrete velocity space in the simulation is proposed and validated through
numerical experiments. Compared with the original UGKS method, the current
adaptive scheme frees the memory requirement in near-equilibrium flow regimes
and speeds up the computation, but provides equivalent physical solutions. Due
to the use of particle distribution function in the whole flow field, the
AUGKS avoids defining an interface condition for domain decomposition in order
to connect fluid and kinetic representations. In other words, no buffer zone
is needed in the AUGKS.

This paper is organized as following. Section 2 is a brief introduction of the
fundamental kinetic theory. Section 3 presents the numerical implementation of
the adaptive unified gas kinetic scheme and proposes a switching criterion of
the velocity space transformation. Section 4 includes numerical examples to
demonstrate the performance of the current scheme. The last section is the
conclusion.

##  2 Gas kinetic modeling

The gas kinetic theory describes the time-space evolution of particle
distribution function $f(\mathbf{x},\mathbf{u},t)$, where
$\mathbf{x}\in\mathcal{R}^{3}$ is space variable and
$\mathbf{u}\in\mathcal{R}^{3}$ is particle velocity. In the absence of
external force field, the Boltzmann equation of a monatomic dilute gas writes,

$f_{t}+\mathbf{u}\cdot\nabla_{\mathbf{x}}f=Q(f,f)=\int_{\mathcal{R}^{3}}\int_{\mathcal{S}^{2}}\left[f(\mathbf{u}^{\prime})f(\mathbf{u}_{1}^{\prime})-f(\mathbf{u})f(\mathbf{u}_{1})\right]B(\cos\theta,g)d\mathbf{\Omega}d\mathbf{u}_{1},$
(1)

where $\mathbf{u},\mathbf{u_{1}}$ are the pre-collision velocities of two
classes of molecules, and $\mathbf{u}^{\prime},\mathbf{u_{1}}^{\prime}$ are
the corresponding post-collision velocities. The collision kernel
$B(\cos\theta,g)$ measures the strength of collisions in different directions,
where $\theta$ is the deflection angle and
$g=|\mathbf{g}|=|\mathbf{u}-\mathbf{u_{1}}|$ is the magnitude of relative pre-
collision velocity. The $\mathbf{\Omega}$ is the unit vector along the
relative post-collision velocity
$\mathbf{u}^{\prime}-\mathbf{u_{1}}^{\prime}$, and the deflection angle
$\theta$ satisfies $\cos\theta=\mathbf{\Omega}\cdot\mathbf{g}/g$. The
conservation of momentum and energy leads the following relations,

$\displaystyle\mathbf{u}^{\prime}=\frac{\mathbf{u}+\mathbf{u_{1}}}{2}+\frac{|\mathbf{u}-\mathbf{u_{1}}|}{2}\mathbf{\Omega}=\mathbf{u}+\frac{g\mathbf{\Omega}-\mathbf{g}}{2},$
(2)
$\displaystyle\mathbf{u_{1}}^{\prime}=\frac{\mathbf{u}+\mathbf{u_{1}}}{2}-\frac{|\mathbf{u}-\mathbf{u_{1}}|}{2}\mathbf{\Omega}=\mathbf{u}_{1}-\frac{g\mathbf{\Omega}-\mathbf{g}}{2}.$

Due to the complicated fivefold integration in the Boltzmann collision
operator, some simplified kinetic models have been constructed, such as the
Shakhov [18]. In this model, the Boltzmann collision operator $Q(f,f)$ is
replaced with a relaxation operator $S(f)$, which writes,

$\displaystyle
f_{t}+\mathbf{u}\cdot\nabla_{\mathbf{x}}f=S(f)=\frac{f^{+}-f}{\tau},$ (3)
$\displaystyle
f^{+}=\rho\left(\frac{\lambda}{\pi}\right)^{\frac{3}{2}}e^{-\lambda(\mathbf{u}-\mathbf{U})^{2}}\left[1+(1-{Pr})\mathbf{c}\cdot\mathbf{q}\left(\frac{\mathbf{c}^{2}}{RT}-5\right)/(5pRT)\right],$

where $\tau=\mu/p$ is the collision time. The macroscopic density, velocity,
temperature and heat flux are marked with $\rho,\mathbf{U},T,\mathbf{q}$. The
$\mathbf{c}=\mathbf{u}-\mathbf{U}$ is particle peculiar velocity, $Pr$ is the
Prandtl number, $R$ is the gas constant, and $\lambda={\rho}/{(2p)}$. In this
paper, the numerical simulations will be conducted by either the full
Boltzmann or the Shakhov collision terms.

The macroscopic conservative flow variables are related with the moments of
particle distribution function via

$\textbf{W}=\left(\begin{matrix}\rho\\\ \rho\mathbf{U}\\\ \rho
E\end{matrix}\right)=\int f\psi d\mathbf{u},$

and the collision terms satisfy the compatibility condition,

$\int Q(f,f_{1})\psi d\mathbf{u}=\int S(f)\psi d\mathbf{u}=0,$

where $\psi=\left(1,\mathbf{u},\frac{1}{2}\mathbf{u}^{2}\right)^{T}$ is a
vector of collision invariants. Here we rewrite the collision terms $Q(f,f)$
and $S(f)$ into a general form $Q(f)$.

With a local constant collision time $\tau$, the integral solution of Eq.(3)
can be constructed along the characteristic line,

$f(\mathbf{x},t,\mathbf{u})=\frac{1}{\tau}\int_{t^{0}}^{t}f^{+}(\mathbf{x}^{\prime},t^{\prime},\mathbf{u})e^{-(t-t^{\prime})/\tau}dt^{\prime}+e^{-t/\tau}f_{0}(\mathbf{x}^{0},0,\mathbf{u}),$
(4)

where $\mathbf{x}^{\prime}=-\mathbf{u}(t-t^{\prime})$ is the particle
trajectory, and $f_{0}$ is the gas distribution function at the initial time
step $t=t^{0}$. Based on the above evolving solution, the corresponding
discretized gas dynamic equations on the cell size and time step scales can be
constructed in the gas kinetic scheme.

##  3 Adaptive unified gas kinetic scheme

In this section, we will present the principle and numerical implementation of
the velocity-space adaptive unified gas kinetic scheme (AUGKS). The original
gas kinetic scheme with continuous and discrete particle velocity space will
be introduced first. The detailed coupling of continuum and kinetic treatments
and the switching criterion for velocity space transformation will be
discussed. For simplicity, the following introduction is based on two-
dimensional case, while the extension to three dimension is straightforward.

###  3.1 Unified gas kinetic scheme with discrete velocity space

With the notation of cell averaged quantities in the control volume,

$\displaystyle\mathbf{W}_{x_{i},y_{j},t^{n}}=\mathbf{W}_{i,j}^{n}=\frac{1}{\Omega_{i,j}(\mathbf{x})}\int_{\Omega_{i,j}}f(x,y,t^{n})dxdy,$
$\displaystyle
f_{x_{i},y_{j},t^{n},u_{l},v_{m}}=f_{i,j,l,m}^{n}=\frac{1}{\Omega_{i,j}(\mathbf{x})\Omega_{l,m}(\mathbf{u})}\int_{\Omega_{i,j}}\int_{\Omega_{l,m}}f(x,y,t^{n},u,v)dxdydudv,$

where $\Omega_{i,j}$ and $\Omega_{l,m}$ are the cell area in the discretized
physical and velocity space, the macroscopic conservative variables and the
particle distribution function are updated in the UGKS,

$\textbf{W}_{i,j}^{n+1}=\textbf{W}_{i,j}^{n}+\frac{1}{\Omega_{i,j}}\int_{t^{n}}^{t^{n+1}}\sum_{r}{\mathbf{F}}_{r}\Delta
L_{r}dt,$ (5) $\displaystyle f_{i,j,l,m}^{n+1}=$ $\displaystyle
f_{i,j,l,m}^{n}+\frac{1}{\Omega_{i,j}}\int_{t^{n}}^{t^{n+1}}\sum_{r}u_{r}f_{r}\Delta
L_{r}dt+\int_{t^{n}}^{t^{n+1}}\int_{\Omega_{i,j}}Q(f)dt,$ (6)

where $\textbf{F}_{r}$ is the flux of conservative variables, $f_{r}$ is the
time-dependent gas distribution function at cell interface and $\Delta L_{r}$
is the cell interface length.

In the UGKS, the flux function $\mathbf{F}_{r}$ is evaluated through the
particle distribution function at the cell interface, which is constructed
from the evolving solution of the Shakhov equation. If the interface location
is simplified with $x_{i+1/2}=0$, $y_{j}=0$, and $t^{n}=0$, with a local
constant collision time $\tau$, the integral solution in Eq.(4) can be written
as,

$\displaystyle f(0,0,t,u_{l},v_{m},\xi)=$
$\displaystyle\frac{1}{\tau}\int_{0}^{t}f^{+}(x^{\prime},y^{\prime},t^{\prime},u_{l},v_{m},\xi)e^{-(t-t^{\prime})/\tau}dt^{\prime}$
(7) $\displaystyle+e^{-t/\tau}f_{0}(x^{0},y^{0},0,u_{l},v_{m},\xi),$

where $x^{\prime}=-u_{l}(t-t^{\prime})$ and $y^{\prime}=v_{m}(t-t^{\prime})$
are the particle trajectories, and $x^{0},y^{0}$ are the initial locations for
the particle which passes through the cell interface at time $t$. Here $f_{0}$
is the gas distribution function at the beginning of $n$-th time step. The
internal degree of freedom $\xi$ denotes the random motion in $z$ direction.
The above scale-dependent evolving solution provides a multiple scale
evolution model, where the contributions from both equilibrium hydrodynamic
and non-equilibrium kinetic flow physics are included uniformly.

In the detailed numerical scheme, to the second order accuracy, the initial
gas distribution function $f_{0}$ is reconstructed as

$f_{0}(x,y,0,u_{l},v_{m},\xi)=\left\\{\begin{aligned}
&f_{i+1/2,j,l,m}^{L}+\sigma_{i,j,l,m}x+\theta_{i,j,l,m}y,\quad x\leq 0,\\\
&f_{i+1/2,j,l,m}^{R}+\sigma_{i+1,j,l,m}x+\theta_{i+1,j,l,m}y,\quad
x>0,\end{aligned}\right.$ (8)

where $f_{i+1/2,j,l,m}^{L}$ and $f_{i+1/2,j,l,m}^{R}$ are the reconstructed
initial distribution functions at the left and right hand sides of a cell
interface, and $\sigma$ and $\theta$ are the slopes of distribution function
along $x$ and $y$ directions. In addition, the equilibrium distribution
function around a cell interface is constructed as

$f^{+}=f^{+}_{0}\left[1+(1-H[x])\bar{a}^{L}x+H[x]\bar{a}^{R}x+\bar{b}y+\bar{A}t\right],$
(9)

where $f^{+}_{0}$ is the equilibrium distribution at $(x=0,y=0,t=0)$. The
coefficients $(\bar{a}^{L,R},\bar{b},\bar{A})$ can be evaluated from the
spatial distribution of conservative variables on both sides of the cell
interface and the compatibility condition [19]. After all coefficients are
determined, the time dependent interface distribution function becomes

$\displaystyle f(0,0,t,u_{l},v_{m},\xi)=$
$\displaystyle\left(1-e^{-t/\tau}\right)f^{+}_{0}$ (10)
$\displaystyle+\left(\tau(-1+e^{-t/\tau})+te^{-t/\tau}\right)\bar{a}^{L,R}u_{l}f^{+}_{0}$
$\displaystyle+\left(\tau(-1+e^{-t/\tau})+te^{-t/\tau}\right)\bar{b}v_{m}f^{+}_{0}$
$\displaystyle+\tau\left(t/\tau-1+e^{-t/\tau}\right)\bar{A}f^{+}_{0}$
$\displaystyle+e^{-t/\tau}\left[\left(f_{i+1/2,j,l,m}^{L}-u_{l}t\sigma_{i,j,l,m}-v_{m}t\theta_{i,j,l,m}\right)H\left[u_{l}\right]\right.$
$\displaystyle\left.+\left(f_{i+1/2,j,l,m}^{R}-u_{l}t\sigma_{i+1,j,l,m}-v_{m}t\theta_{i+1,j,l,m}\right)(1-H\left[u_{l}\right])\right]$
$\displaystyle=$
$\displaystyle\widetilde{f}^{+}_{i+1/2,j,k,l}+\widetilde{f}_{i+1/2,j,k,l},$

where $\widetilde{f}^{+}_{i+1/2,j,k,l}$ is related to equilibrium state
integration and $\widetilde{f}_{i+1/2,j,k,l}$ is related to the initial
distribution. With the variation of the ratio between evolving time $t$ (i.e.,
the time step in the computation) and collision time $\tau$, the interface
distribution function above provides self-conditioned multiple scale flow
physics across different flow regimes. After the interface distribution
function is determined, the corresponding fluxes of conservative variables are
evaluated through

${\textbf{F}}_{r}=\int_{\Omega_{l,m}}u_{l}f(0,0,t,u_{l},v_{m},\xi)\psi
dudvd\xi.$

Inside each control volume, the collision term $Q(f)$ is to be determined for
the update of particle distribution function in Eq.(6). In the unified scheme,
the numerical treatments of $Q(f)$ for the exact Boltzmann collision term or
the Shakhov model are given in [20, 21, 22, 23].

###  3.2 Near-continuum gas kinetic scheme with continuous velocity space

In continuum flow with intensive intermolecular collisions, the particle
distribution function is close to local equilibrium state, and thus the
Navier-Stokes equations are valid to describe macroscopic fluid motion. In
this case, it is straightforward to apply the Chapman-Enskog expansion to
construct the corresponding distribution function at every time step, and thus
only macroscopic conservative variables need to be stored and updated.
Following this procedure, in the gas kinetic scheme with continuous particle
velocity space [24], the gas distribution function $f_{0}$ in Eq.(7) at the
beginning of each time step is expanded as following,

$f_{0}=\left\\{\begin{aligned}
&f_{0}^{+(l)}\left(1+a^{l}x+b^{l}x-\tau\left(a^{l}u+b^{l}v+A^{l}\right)\right),\quad
x\leq 0\\\
&f_{0}^{+(r)}\left(1+a^{r}x+b^{r}x-\tau\left(a^{r}u+b^{r}v+A^{r}\right)\right),\quad
x>0\end{aligned}\right.$ (11)

where $f_{0}^{+(l)}$ and $f_{0}^{+(r)}$ are the corresponding equilibrium
state at left and right sides of the interface, and the coefficients
$(a^{l,r},b^{l,r},A^{l,r})$ are related with spatial distribution of
macroscopic variables. The equilibrium distribution function $f^{+}$ around a
cell interface is constructed in the same way in Eq.(9), and the corresponding
interface distribution function becomes,

$\displaystyle f(0,0,t,u,v,\xi)=$
$\displaystyle\left(1-e^{-t/\tau}\right)f_{0}^{+}$ (12)
$\displaystyle+\left(\tau(-1+e^{-t/\tau})+te^{-t/\tau}\right)(\bar{a}^{L,R}u+\bar{b}v)f^{+}_{0}$
$\displaystyle+\tau\left(t/\tau-1+e^{-t/\tau}\right)\bar{A}f^{+}_{0}$
$\displaystyle+e^{-t/\tau}\left\\{[1-(t+\tau)(ua^{l}+vb^{l})]H[u]f_{0}^{+(l)}\right.$
$\displaystyle\left.+[1-(t+\tau)(ua^{r}+vb^{r})][1-H[u]]f_{0}^{+(r)}\right\\}$
$\displaystyle+e^{-t/\tau}\left[-\tau A^{l}H[u]f_{0}^{+(l)}-\tau
A^{r}(1-H[u])f_{0}^{+(r)}\right].$

The interface distribution function here is in the continuous form of particle
velocity $(u,v)$, and the fluxes for macroscopic variables can be obtained
through the analytical integration of Gaussian distribution.

###  3.3 Adaptive unified gas kinetic scheme

In a multiscale flow problem, to overcome the computational deficiency and
memory burden from a large amount of discretized velocity points, it is
feasible to combine both continuum and rarefied flow solvers into a single
framework with dynamic continuous-discrete velocity transformation. As shown
in Fig.1, in near-equilibrium flow regions, the particle distribution function
(PDF) is formulated with a continuous velocity space based on the Chapman-
Enskog expansion, and only macroscopic flow variables are updated. In non-
equilibrium regions, the AUGKS switches to a discretized velocity space to
follow the evolution of particle distribution function. The continuous and
discrete velocity space are connected with an adaptation interface, at which
the continuous solution of distribution function is sorted into discretized
velocity space.

In the detailed numerical scheme, the macroscopic conservative variables are
updated in Eq.(5), while in non-equilibrium regions the particle distribution
function is updated in Eq.(6). Near the adaptation interface, at every time
step $t^{n}$, if there is no recorded discretized distribution function at
$t^{n-1}$ in the current ”non-equilibrium” cell $(i,j)$, a local discretized
velocity mesh will be generated, where the particle distribution function at
velocity point $(l,m)$ is given by the discrete Chapman-Enskog expansion,

$f_{i,j,l,m}=f^{+}_{0(i,j,l,m)}[1-\tau(au_{l}+bv_{m}+A)],$ (13)

where the coefficients $a,b,A$ are related with the spatial distribution of
macroscopic variables, and can be determined in the same way in Sec. 3.2. In
the current scheme, the velocity mesh is generated within

$u\in[U-3\sqrt{RT},U+3\sqrt{RT}],\ v\in[V-3\sqrt{RT},V+3\sqrt{RT}],$ (14)

where $(U,V)$ is macroscopic flow velocity, $T$ is temperature and $R$ is the
gas constant.

To update the discretized distribution function in the adjacent cell next to
the adaptation interface, the interface distribution function from the
continuous GKS solution in Eq.(12) is rewritten into the following discrete
form,

$\displaystyle f(0,0,t,u_{l},v_{m},\xi)=$
$\displaystyle\left(1-e^{-t/\tau}\right)f_{0}^{+}$ (15)
$\displaystyle+\left(\tau(-1+e^{-t/\tau})+te^{-t/\tau}\right)(\bar{a}^{L,R}u_{l}+\bar{b}v_{m})f^{+}_{0}$
$\displaystyle+\tau\left(t/\tau-1+e^{-t/\tau}\right)\bar{A}f^{+}_{0}$
$\displaystyle+e^{-t/\tau}\left\\{[1-(t+\tau)(u_{l}a^{l}+v_{m}b^{l})]H[u_{l}]f_{0}^{+(l)}\right.$
$\displaystyle\left.+[1-(t+\tau)(u_{l}a^{r}+v_{m}b^{r})][1-H[u_{l}]]f_{0}^{+(r)}\right\\}$
$\displaystyle+e^{-t/\tau}\left[-\tau A^{l}H[u_{l}]f_{0}^{+(l)}-\tau
A^{r}(1-H[u_{l}])f_{0}^{+(r)}\right],$

and then it can be used to determine the fluxes of macroscopic variables and
particle distribution function. In this way, the fluxes at adaptation
interface are fully determined and can be used to update the macroscopic
variables in Eq.(5) and particle distribution function in Eq.(6).

In the current scheme, the time step is determined by the CFL condition,

$\Delta t=\mathrm{CFL}\frac{\Delta x\Delta y}{u_{max}\Delta y+v_{max}\Delta
x},$ (16)

where CFL is the CFL number, and $(u_{max},v_{max})$ is the largest particle
velocity in $x$ and $y$ directions.

$\underbrace{\hskip 196.32384pt}$$\underbrace{\hskip 196.32384pt}$velocity
space adaptationnear-equilibrium regionnon-equilibrium regioncontinuous
velocitydiscrete velocity Figure 1: Schematic of the adaptive scheme for
multiscale flow.

###  3.4 Switching criterion of velocity space

The accuracy and efficiency of the current adaptive scheme are based on a
proper choice of location of velocity space adaptation. The transition of
discrete to continuous velocity space must be located in the region where the
Navier-Stokes solutions provided by the GKS with a continuous velocity space
are still valid. For this issue, many empirical parameters for the breakdown
of continuum description have been proposed. Bird [4] proposed a parameter
$P=D(ln\rho)/Dt/\nu$ for the DSMC simulation of expansion flows, where $\rho$
is gas density and $\nu$ is collision frequency, and the breakdown value of
$P$ for translational equilibrium is $0.05$. Boyd et al. [25] extended the
above concept to a gradient-length-local Knudsen number $Kn_{GLL}=\ell|\nabla
Q|/{Q}$, where $\ell$ is the local mean free path and $Q$ is the macroscopic
flow quantity of interest, with the critical value $0.05$. Considering the
terms in the Chapman-Enskog distribution function, Garcia et al. [30] proposed
a breakdown parameter based on dimensionless stress and heat flux
$B=max(|\tau^{*}|,|q^{*}|)$, with the switching criterion of $0.1$.

Since the particle distribution function takes the Chapman-Enskog expansion
form in the evolution process of the continuous GKS solver, here we propose an
alternative switching criterion of particle velocity space directly from the
Chapman-Enskog expansion. For brevity, one-dimensional case will be used for
illustration. When there is no discontinuity inside the flow field, the
particle distribution function in the Chapman-Enskog form writes,

$f=f_{0}^{+}[1-\tau(au+A)],$ (17)

where the coefficients $a$ and $A$ are the space and time variations of a
distribution function, which can be expanded based on the collision
invariants,

$\displaystyle a=a_{1}+a_{2}u+a_{3}\frac{1}{2}u^{2}=a_{\alpha}\psi_{\alpha},$
$\displaystyle A=A_{1}+A_{2}u+A_{3}\frac{1}{2}u^{2}=A_{\alpha}\psi_{\alpha}.$

In the current one-dimensional case, the spatial slope $a$ is determined from
the following relation,

$\frac{\partial\mathbf{W}}{\partial x}=\int af^{+}_{0}\psi
d\mathbf{u}=M\mathbf{a},$

where $\psi$ is a vector of collision invariants, $M_{\alpha,\beta}=\int
f^{+}_{0}\psi_{\alpha}\alpha\psi_{\alpha}\beta du$ and
$\mathbf{a}=(a_{1},a_{2},a_{3})^{T}$. The solution $\mathbf{a}$ can be written
as,

$\displaystyle a_{3}=4\frac{\lambda^{2}}{(K+1)\rho}\left[2\frac{\partial(\rho
E)}{\partial x}-2U\frac{\partial(\rho U)}{\partial
x}+\frac{\partial\rho}{\partial
x}\left(U^{2}-\frac{K+1}{2\lambda}\right)\right],$ $\displaystyle
a_{2}=2\frac{\lambda}{\rho}\left[2\frac{\partial(\rho U)}{\partial
x}-U\frac{\partial\rho}{\partial x}\right]-Ua_{3},$ $\displaystyle
a_{1}=\frac{1}{\rho}\frac{\partial\rho}{\partial
x}-Ua_{2}-\frac{1}{2}\left[U^{2}+\frac{K+1}{2\lambda}\right]a_{3}.$

The time derivative $A$ is related to the temporal variation of conservative
flow variables respectively, and can be evaluated from a successive procedure
in the Chapman-Enskog expansion with the help of the Euler equations,

$\frac{\partial\mathbf{W}}{\partial t}=\int Af^{+}_{0}\psi d\mathbf{u}=-\int
u\frac{\partial f_{0}^{+}}{\partial x}\psi d\mathbf{u}.$

Generally, the Navier-Stokes solutions can be applied when the Chapman-Enskog
expansion is a proper approximation of the distribution function in near-
equilibrium regime. Therefore, based on the dimensionless collision time and
space-time variations, a switching criterion for the velocity space
transformation can be defined as

$B=\tau\mathrm{max}(|{\hat{a}}|,|{\hat{A}}|),$ (18)

where dimensionless variables are defined as,

$\hat{\tau}=\frac{\tau(2RT_{0})^{1/2}}{L_{0}},{\hat{a}}=a{L_{0}},{\hat{A}}=\frac{AL_{0}}{(2RT_{0})^{1/2}},$
(19)

where $R$ is gas constant, $L_{0}$ and $T_{0}$ are reference length and
temperature. The current switching criterion $B$ for particle velocity space
will be tested in numerical experiments.

###  3.5 Summary

The numerical algorithm of the adaptive unified gas kinetic scheme is as
following. In the AUGKS, we follow the evolution of both both conservative
variables and particle distribution function. In near-equilibrium flow
regions, the particle distribution function is formulated by the Chapman-
Enskog expansion in a continuous velocity space, and macroscopic flow
variables are updated in Eq.(5). For non-equilibrium flows, the particle
distribution function is updated explicitly in Eq.(6). The scale-dependent
flux function is determined by the particle distribution function at the
interface, which comes from the integral solution of kinetic model equation in
Eq.(4). In each time step, the domain of continuous and discrete velocity
space is specified by Eq.(18), and the corresponding interface fluxes are
provided by Eq.(10), Eq.(15) and Eq.(12). The detailed numerical procedures
for AUGKS are given in Fig. 2.

StartCalculate time step by Eq.(16)Evaluate velocity adaptation by
Eq.(18)Continuous velocityDiscrete velocityCalculate flux by Eq.(12)Calculate
flux by Eq.(15)Calculate flux by Eq.(10)destroy recorded PDFConstruct
unrecorded PDF by Eq.(13)Update conservative variables by Eq.(5)Calculate
equilibrium distribution $f^{+(n+1)}$ and collision time $\tau^{n+1}$Update
PDF by Eq.(6)Update conservative variables by Eq.(5)end Figure 2: Numerical
algorithm of AUGKS.

##  4 Numerical experiments

In this section, we are going to present some numerical experiments to test
the performance of the current AUGKS. In order to demonstrate multiscale
performance of the algorithm, simulations from continuum Euler and Navier-
Stokes to free molecule flow are presented. The following dimensionless flow
variables are used in the calculations,

$\displaystyle\hat{x}=\frac{x}{L_{0}},\hat{y}=\frac{y}{L_{0}},\hat{\rho}=\frac{\rho}{\rho_{0}},\hat{T}=\frac{T}{T_{0}},\hat{u}_{i}=\frac{u_{i}}{(2RT_{0})^{1/2}},$
$\displaystyle\hat{U}_{i}=\frac{U_{i}}{(2RT_{0})^{1/2}},\hat{f}=\frac{f}{\rho_{0}(2RT_{0})^{3/2}},\hat{P}_{ij}=\frac{P_{ij}}{\rho_{0}(2RT_{0})},\hat{q}_{i}=\frac{q_{i}}{\rho_{0}(2RT_{0})^{3/2}},$

where $u_{i}$ is the particle velocity, $U_{i}$ is the macroscopic flow
velocity, $P_{ij}$ is the stress tensor, $q_{i}$ is the heat flux. The
subscript zero represents the reference state. For simplicity, the hat
notation to denote dimensionless variables will be removed henceforth. Argon
gas is used in the simulation with the variable hard sphere (VHS) molecule
model, and the dynamic viscosity is related with the Knudsen number in the
reference state via

$\mu_{ref}=\frac{5(\alpha+1)(\alpha+2)\sqrt{\pi}}{4\alpha(5-2\omega)(7-2\omega)}Kn_{ref}.$

In this simulation, we choose $\alpha=1.0$ and $\omega=0.5$ to recover a hard
sphere gas, and the viscosity varies with temperature through

$\mu=\mu_{ref}\left(\frac{T}{T_{ref}}\right)^{\theta},$

where $T_{ref}$ is the reference temperature and $\theta=0.81$ is the index of
viscosity coefficient.

###  4.1 Shock tube problem

The first case is the Sod shock tube problem. The flow domain $x\in[0,1]$ is
divided into 100 uniform cells. The initial condition is set as

$\displaystyle\rho=1.0,U=0.0,p=1.0,\ x\leq 0.5,$
$\displaystyle\rho=0.125,U=0.0,p=0.1,\ x>0.5.$

The simulations are performed with reference Knudsen numbers varying from
$Kn=0.0001$ to $Kn=1.0$, corresponding to different flow regimes. The current
criterion value for velocity space transformation is set as $B=0.0001$. The
velocity space is discretized into $80$ uniform points for the update of
particle distribution function. In the AUGKS and UGKS, the full Boltzmann
collision operator is solved here by the fast spectral method [21]. The
numerical solutions at $t=0.2$ are presented in Fig. 3, 4, 5 and 6. The
reference solution of continuum flow is calculated by the continuous GKS
solver with 1000 cells, and the free molecular flow solution is derived from
the collisionless Boltzmann equation.

In the simulation, the region with initial homogeneous spatial distribution of
flow variables is calculated with a continuous velocity space, except the
central discontinuity simulated with a discretized velocity. As time evolves,
the non-equilibrium region enlarges along with the use of discrete velocity
space. As presented in Fig. 7a, at $Kn=0.0001$ and $t=0.2$, the flow domain is
divided into some subzones, where the non-equilibrium particle distribution
function inside rarefaction wave, contact discontinuity, and shock wave is
fully resolved with the discretized velocity space, while in the rest near-
equilibrium regions the Chapman-Enskog expansion is adopted over a continuous
velocity space. The solutions of AUGKS at $Kn=0.0001$ and $t=0.2$ are
presented in Fig. 3a, 3b, 3c, which match the benchmark continuum and UGKS
solutions accurately. As the Knudsen number gets to $Kn=0.001$ at $t=0.2$,
near-equilibrium region confines to a small part near the left tube boundary,
where the distribution function has a continuous velocity space, which is
shown in Fig. 7b. With increasing rarefaction effect, the distributions of
flow variables deviate from the NS solutions gradually and tend to
collisionless Boltzmann solutions. As the reference Knudsen number gets to
$Kn=0.01$, the Navier-Stokes solutions lose its validity quickly from the
initial condition, and the non-equilibrium region occupies the whole tube at
$t=0.2$. The numerical solution approaches to the collisionless Boltzmann
solution at $Kn=1.0$, as shown in Fig. 5 and 6.

This test case illustrates the capacity of AUGKS to simulate flow in different
regimes. The asymptotic preserving (AP) property is confirmed in the two
limiting solutions. With increasing reference Knudsen number, there is a
smooth transition from the Euler solution of the Riemann problem to
collisionless Boltzmann solution. Table 1 presents the CPU time cost from the
current adaptive scheme and original UGKS method. As is shown, the AUGKS is
about $3.63$ times faster than UGKS at $Kn=0.0001$. Since the computational
cost is proportional to the mesh points in the velocity space, it is expected
that the computational efficiency is closely related to the size of non-
equilibrium region. When the rarefaction increases, the CPU time of AUGKS
increases correspondingly, while it is still more efficient than the original
UGKS.

| AUGKS | UGKS
---|---|---
Kn=0.0001 | 2042.98 | 7421.70
Kn=0.001 | 3537.73 | 7430.82
Kn=0.01 | 4692.52 | 7547.99
Kn=0.1 | 5694.46 | 7275.40
Table 1: CPU time cost in the Sod case using AUGKS and original UGKS at
different Knudsen numbers.

(a) Density (b) Velocity (c) Temperature

Figure 3: Sod shock tube at $t=0.2$ with reference Knudsen number $Kn=0.0001$.

(a) Density (b) Velocity (c) Temperature

Figure 4: Sod shock tube at $t=0.2$ with reference Knudsen number $Kn=0.001$.

(a) Density (b) Velocity (c) Temperature

Figure 5: Sod shock tube at $t=0.2$ with reference Knudsen number $Kn=0.01$.

(a) Density (b) Velocity (c) Temperature

Figure 6: Sod shock tube at $t=0.2$ with reference Knudsen number $Kn=1.0$.

(a) $Kn=0.0001$ (b) $Kn=0.001$

Figure 7: Velocity space adaptation inside the shock tube at $t=0.2$.

###  4.2 Rayleigh flow

A Rayleigh flow forms over a plate which suddenly acquires a constant parallel
velocity and temperature. In this test case, we follow the setup by Sun [26].
As shown in Fig. 8, the argon gas is at rest and has a unit temperature
initially. When $t>0$, the plate suddenly moves with a constant velocity
$U_{w}=0.0296$ and temperature $T_{w}=1.36$. The momentum and energy are
transported into the flow field through a shearing effect in the unsteady
process. A physical domain $y\in[0,1]$ with 100 uniform cells are set up for
the simulation, and the 32 uniform points are used in the velocity space where
the particle distribution function is updated directy. In this case, the full
Boltzmann collision operator in the AUGKS and UGKS is calculated by the fast
spectrum method [21]. The current switching criterion of velocity space is set
as $B=0.0001$.

Numerical simulations are performed with a series of reference Knudsen number,
and solutions at same output times are plotted in Fig. 9, 10, 11 and 12.
Besides AUGKS solutions, the UGKS and DSMC results are also provided as
benchmarks. With Maxwell’s fully accommodation boundary condition, Bird [4]
proposed an analytical solution from the collisionless Boltzmann equation when
the time is much less than the reference mean collision time
$\tau_{0}=\ell_{0}/v_{0}$, where $\ell_{0}$ is particle mean free path and
$v_{0}$ is the mean molecular speed. The analytical collisionless solution is
also plotted in figures.

As presented in Fig. 9, 10, 11 and 12, for the case at $Kn=2.66$ and
$t=0.1\tau_{0}$, the AUGKS recovers exact collisionless Botlzmann solution. In
the transition regime $Kn=0.266$ and $Kn=0.0266$ at $t=\tau_{0}$ and
$t=10\tau_{0}$, the numerical solutions deviate from collisionless solutions
gradually due to increasing intermolecular collisions. At $t=100\tau_{0}$ and
$Kn=0.00266$ corresponding to a near-continuum regime, the current adaptive
scheme recovers the Navier-Stokes solutions with intensive intermolecular
collisions. As plotted, in all cases the AUGKS solutions match up well with
the benchmark solutions from DSMC and UGKS. It is worth mentioning that in
comparison with DSMC method, the current Boltzmann-equation-based adaptive
unified scheme has no statistical scattering, which is beneficial in low speed
simulations.

Fig. 13 presents the velocity space adaptation inside the flow domain at the
output time. In the case with $Kn=0.00266$, in the near-wall region with large
slope of macroscopic variables, the AUGKS uses a discrete velocity space,
while a continuous velocity space is used in the outer domain. As the Knudsen
number increases, the enhanced dimensionless viscosity and heat conductivity
lead to a large non-equilibrium region. As a result, the non-equilibrium
region enlarges faster. For the case $Kn=0.266$, in all flow region the
distribution function deviates from the Chapman-Enskog solution and its
evolution must be followed with a discretized velocity space. Table 2 presents
the computational cost of the AUGKS and UGKS. When $Kn=0.00266$, the AUGKS is
$7.45$ times faster than the original UGKS. When the Knudsen number increases,
the enhanced non-equilibrium effects increase the computational cost in AUGKS.
From the current numerical experiments, it is clear that the AUGKS provides a
self-adjusted algorithm from continuum to rarefied flow simulation with the
consideration of both accuracy and efficiency.

$T_{w}=1.36,U_{w}=0.0296$$U_{w}$$\rm Outer\ boundary$$\rm Argon\
gas$$T_{0}=1,U_{0}=V_{0}=0$$x$$y$ Figure 8: Schematic of Rayleigh problem.

(a) Density and Temperature (b) Velocity

Figure 9: Rayleigh flow at $t=0.1\tau_{0}$ with reference Knudsen number
$Kn=2.66$.

(a) Density and Temperature (b) Velocity

Figure 10: Rayleigh flow at $t=\tau_{0}$ with reference Knudsen number
$Kn=0.266$.

(a) Density and Temperature (b) Velocity

Figure 11: Rayleigh flow at $t=10\tau_{0}$ with reference Knudsen number
$Kn=0.0266$.

(a) Density and Temperature (b) Velocity

Figure 12: Rayleigh flow at $t=100\tau_{0}$ with reference Knudsen number $Kn=0.00266$. Figure 13: Velocity space adaptation in the Rayleigh flow at the output instant with different reference Knudsen numbers. | AUGKS | UGKS
---|---|---
Kn=0.00266 | 1003.59 | 7477.10
Kn=0.0266 | 2223.62 | 7460.21
Kn=0.266 | 3751.75. | 7438.06
Kn=2.66 | 4637.19 | 7457.87
Table 2: CPU time cost in the Rayleigh problem using the AUGKS and UGKS at
different Knudsen numbers.

###  4.3 Nozzle flow

The nozzle flow connecting different flow regimes is an ideal case to test the
capacity of AUGKS in capturing multiple scale flow dynamics. The schematic of
the nozzle problem is presented in Fig. 14. The argon gas is enclosed in a
rectangular box $x\in[0,2.2],y\in[-0.5,0.5]$. The velocity space is
discretized into $28\times 28$ uniform points for the update of particle
distribution function. The switching criterion of velocity space in this case
is set as $B=0.0005$. In this case, the collision term in AUGKS and UGKS is
the Shakhov model. The computational domain is divided into two parts which
are connected through a nozzle. The gas density in the left is $100$ times
higher than taht in the right part, and the initial Knudsen numbers are
$Kn_{L}=0.0001$ and $Kn_{R}=0.01$ respectively. Except that, the initial gas
is at rest and has the same temperature inside two subdomains, same as the
cavity wall. All the walls are set up as Maxwell’s diffusive boundary. The
nozzle has two entrances with cross sections $\Delta y_{L}=0.13$ and $\Delta
y_{R}=0.33$ along a length $\Delta x=0.14$, from which a jet flow will be
formed. The simulation is performed till $t=50\tau_{0}$, where
$\tau_{0}=\ell_{0}/v_{0}$ is the mean collision time of initial argon gas in
the right domain.

Fig. 15 and 16 presents the solution contours of $U-$velocity and temperature
at two times $t=5\tau_{0}$ and $t=20\tau_{0}$. The upper part of color
contours are the results calculated by the AUGKS (flood) and UGKS (lines),
while the lower part is the Navier-Stokes solutions provided by the GKS with a
continuous velocity space. As is shown, the bow shock and expansion cooling
region behind shock are captured by the two methods. However, it is clear that
at $Kn=0.01$, the Navier-Stokes equations lose validity to quantitatively
describe the flow evolution in the right domain, and it is necessary to use
kinetic method to get accurate solutions here. Fig. 17 and 18 presents the
solutions along the horizontal center line of the box at times $t=5\tau_{0}$
and $t=20\tau_{0}$. It is clear that the AUGKS provide equivalent solutions
with the NS ones in the near-equilibrium left region, and with the Boltzmann
solutions in the non-equilibrium right region. This test demonstrates the
multiscale capability of the adaptive method to get the physical solutions in
the corresponding flow regimes.

Fig. 19, 20 and 21 present the components of spatial slope $a$ used in the
velocity space switching criterion, the mean collision time and the
corresponding velocity space adaptation of the flow domain at three times
$t=\tau_{0},\ 5\tau_{0},\ 20\tau_{0}$. As shown, the shock and expansion waves
are the major sources for large flow gradients inside the domain. Accompanying
with the high-density jet into the right domain, the mean collision time
decreases in the jet region. With time increasing, the local flow structure
becomes more complicated, leading to a large non-equilibrium region. As a
result, the Chapman-Enskog expansion fails in the places where the strong non-
equilibrium effects appear, and the discretized velocity space has to be used
in AUGKS. Table 3 presents the computational cost of AUGKS, UGKS and GKS in
this case. With the current setup of physical mesh and velocity space, the
continuous GKS solver is about 25 times faster than the UGKS with discretized
velocity space. The AUGKS is $1.7$ times faster than the original UGKS. At the
initial stage in the simulation, the memory cost in the AUGKS is on the same
order as the GKS, which is about $1/30$ of the UGKS. As flow evolves, the
number of cells associated with discretized velocity space increases, and the
corresponding memory cost becomes larger. At the final time $t=50\tau_{0}$,
there are about 429 out of 914 total cells using continuous velocity space,
and the corresponding memory burden is $53\%$ of the original UGKS.

$T_{w}=1,U_{w}=0$$\rm wall\
boundary$$Kn=0.0001$$Kn=0.01$$T_{0}=1,U_{0}=V_{0}=0$$T_{0}=1,U_{0}=V_{0}=0$$x$$y$
Figure 14: Schematic of Nozzle jet problem.

(a) U-velocity (b) Temperature

Figure 15: Density and Temperature contours at $t=5\tau_{0}$ (upper flood:
AUGKS, upper lines: UGKS, lower: GKS).

(a) U-velocity (b) Temperature

Figure 16: Density and Temperature contours at $t=20\tau_{0}$ (upper flood:
AUGKS, upper lines: UGKS, lower: GKS).

(a) Density (b) U-velocity (c) Temperature

Figure 17: Solutions along the horizontal central line at $t=5\tau_{0}$.

(a) Density (b) U-velocity (c) Temperature

Figure 18: Solutions along the horizontal central line at $t=20\tau_{0}$. | AUGKS | UGKS | GKS
---|---|---|---
CPU time (s) | 526.25 | 898.69 | 35.25
Memory (KB) $(t=0)$ | 3614 | 74828 | 2636
Memory (KB) $(t=50\tau_{0})$ | 35125 | 74836 | 2640
Table 3: CPU time and memory cost in the nozzle flow at $t=50\tau_{0}$.

(a) $a_{1}$ (b) Colliison time (c) Velocity adaptation

Figure 19: Velocity space adaptation in the nozzle flow at $t=\tau_{0}$.

(a) $a_{1}$ (b) Colliison time (c) Velocity adaptation

Figure 20: Velocity space adaptation in the nozzle flow at $t=5\tau_{0}$.

(a) $a_{1}$ (b) Colliison time (c) Velocity adaptation

Figure 21: Velocity space adaptation in the nozzle flow at $t=20\tau_{0}$.

###  4.4 Flow around circular cylinder

In the previous studies, the AUGKS in the simulation of unsteady flows is well
validated. In this case, the flow passing through a circular cylinder is used
to test the performance of current adaptive scheme for steady flow. The
incoming gas has a uniform velocity with Mach number $Ma=5$ and the same
temperature as the cylinder wall. The reference Knudsen number is set up as
$Kn=0.001$ and $Kn=0.01$ relative to cylinder radius, and the corresponding
dynamic viscosity is $\mu_{ref}=7.313\times 10^{-4}$ and
$\mu_{ref}=7.313\times 10^{-4}$. In the calculation, $60$ cells in radial
direction and $100$ cells in circumferential direction are used in physical
domain. The velocity space is discretized into $41\times 41$ velocity points
for the update of particle distribution function. In this case, the collision
term for particle distribution function in AUGKS and UGKS is modeled as the
Shakhov model. The cylinder’s wall is set up as Maxwell’s diffusive boundary.
The switching criterion of particle velocity space is set as $B=0.0005$.

For steady problem, the computational time can be further reduced with the
help of the GKS with a continuous velocity space. A convergent coarse flow
field can be first obtained by the GKS, and then used as the initial state in
the subsquent adaptive method. The method for the computation of steady flow
is the following.

  1. 1.

From initial setup, use the GKS solver in the entire domain and obtain a
convergent flow field.

  2. 2.

Use the calculated macroscopic flow variables as the initial flow condition to
get the particle distribution function with the discretized Chapman-Enskog
form by Eq.(13).

  3. 3.

Adapt the velocity space based on the current switching criterion in Eq.(18).

  4. 4.

Use continuous velocity space in near-equilibrium flow region and discretized
velocity space in non-equilibrium one, and redo the velocity adaptation after
some iterations until a convergent flow field is obtained.

Fig. 22 and Fig. 23 presents the solution contours of $U-$velocity and
temperature calculated by the AUGKS, UGKS and GKS methods respectively around
the cylinder. The upper part of contours is the results of AUGKS (flood) and
UGKS (solid line) solutions, and the lower part is the GKS solutions. As
shown, the bow shock and expansion cooling region behind shock are well
captured by all methods.

Fig. 24, 25, 26, 27 present the solutions along the horizontal center line in
front of and behind the cylinder. At $Kn=0.001$, the cell size and time step
in the computation are much larger than particle mean free path and collision
time. Due to the limited time-space resolution, all the three methods become
shock-capturing schemes, and a sharp shock profile is obtained in front of the
cylinder in Fig. 24. Near the cylinder wall, due to the non-equilibrium gas
dynamics in gas-surface interaction, there is a slight difference in the
solutions provided by UGKS and GKS. At the same time, the gas density reduces
much in the wake region behind cylinder with emerging rarefied effects, and
there is a significant difference between UGKS and GKS solutions in Fig. 25.

When the reference Knudsen number gets to $Kn=0.01$, the increasing particle
mean free path collision leads to a wider shock structure. This non-
equilibrium evolution is provided in the scale-dependent interface solution
used in AUGKS and UGKS. However, the Chapman-Enskog expansion can only provide
incomplete information about this process in continuous GKS solver. As a
result, the GKS presents a narrower shock profile than that in AUGKS and UGKS
in Fig. 26. In the wake region, due to the increasing collision time at
$Kn=0.01$, the results provided by GKS differ significantly from AUGKS and
UGKS solutions in Fig. 27. The GKS with a continuous velocity space fails to
predict physical solutions in these regions, and particle distribution
function is updated explicitly in the AUGKS with a discretized particle
velocity space. It is clear that the current velocity adaptive unified scheme
captures the equivalent physical solutions as the NS in the near-equilibrium
region and the UGKS ones in the non-equilibrium region.

Fig. 28 and 29 presents two components of spatial slope $\mathbf{a}$ used in
the velocity-space switching criterion, the mean collision time and the
adaptation of velocity space. As can be seen, the shock wave and boundary are
two sources for high gradients of flow variables, leading to the failure of
Chapman-Enskog expansion and Navier-Stokes solutions. Behind the cylinder, the
low density wake leads to an increasing collision time, which is shown in Fig.
28c and Fig. 29c. Therefore, a velocity adaptation is determined as shown in
Fig. 28d and Fig. 29d . The incoming flow as well as a small region between
the bow shock and cylinder is computed with continuous velocity space, while
the rest non-equilibrium region are simulated by the UGKS with discrete
velocity space. Due to increasing collision time and non-equilibrium flow
dynamics, the region under discrete particle velocity space enlarges at
$Kn=0.01$ than that at $Kn=0.01$.

Table 4 and 5 presents the computational cost of AUGKS, UGKS and GKS at
$Kn=0.001$ and $Kn=0.01$ respectively. With the current setup of physical mesh
and velocity space, the continuous GKS solver is about 30 times faster than
the UGKS, and the AUGKS is about $3.3$ times faster than the original UGKS in
this steady flow problem. In the convergent steady state, there are about 3196
at $Kn=0.001$ and 1992 at $Kn=0.01$ out of 6000 total cells using continuous
velocity space, and the corresponding memory burden is about $47\%$ and $67\%$
of the original UGKS.

| AUGKS | UGKS | GKS
---|---|---|---
CPU time (s) | 36130.68 | 117371.67 | 2975.07
Memory (KB) | 452508 | 857520 | 14652
Table 4: CPU time and memory cost in the flow around circular cylinder at $Kn=0.001$. | AUGKS | UGKS | GKS
---|---|---|---
CPU time (s) | 22145.10 | 75510.33 | 2536.55
Memory (KB) | 614542 | 856944 | 12636
Table 5: CPU time and memory cost in the flow around circular cylinder at
$Kn=0.01$.

##  5 Conclusion

The gas dynamics has intrinsic multiple scale nature due to the large
variations of density and characteristic length scale of the flow structures.
Based on scale-dependent time evolving solution of the Boltzmann model
equation, a velocity-space adaptive unified gas kinetic scheme has been
developed in this paper for the simulation of multiscale flow transport. The
current adaptive algorithm is based on a dynamic velocity-space
transformation, where the particle velocity space is continuous in near-
equilibrium region and discrete in non-equilibrium one. A switching criterion
for particle velocity space transformation is proposed based on the Chapman-
Enskog expansion and then is validated through numerical experiments. With a
unified framework with the adaptation of particle velocity space only, the
AUGKS needs no buffer zone for the connection between continuum and kinetic
solutions. This compact property leads to an effective method for multiscale
flow simulation with unsteadiness and complex geometries. Compared with
single-velocity-space framework of the original UGKS, the AUGKS is more
efficient and less memory demanding for multiscale flow computations. The
AUGKS provides a useful tool for non-equilibrium flow studies, and it can be
further improved in the future with the combination of implicit and multigrid
techniques [27, 28].

## Acknowledgement

The authors would like to thank Dr. Chang Liu and Dr. Lei Wu for the help on
numerical implementation of the fast spectral method for the Boltzmann
collision term. The current research is supported by Hong Kong research grant
council (16207715, 16211014, 16206617), and National Science Foundation of
China (11772281,91530319).

## References

  * [1] Sydney Chapman and Thomas George Cowling.  The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases.  Cambridge university press, 1970.
  * [2] Hsue-Shen Tsien.  Superaerodynamics, mechanics of rarefied gases.  Journal of the Aeronautical Sciences, 13(12):653–664, 1946.
  * [3] Vasiliĭ Aristov.  Direct methods for solving the Boltzmann equation and study of nonequilibrium flows.
  * [4] Graeme Austin Bird.  Molecular gas dynamics and the direct simulation of gas flows.  1994\.
  * [5] Jean-François Bourgat, Patrick Le Tallec, and Moulay Tidriri.  Coupling boltzmann and navier–stokes equations by friction.  Journal of Computational Physics, 127(2):227–245, 1996.
  * [6] S Tiwari.  Coupling of the boltzmann and euler equations with automatic domain decomposition.  Journal of Computational Physics, 144(2):710–726, 1998.
  * [7] Iain D Boyd and Timothy R Deschenes.  Hybrid particle-continuum numerical methods for aerospace applications.  Technical report, Michigan Univ Ann Arbor Dept of Aerospace Engineering, 2011.
  * [8] HS Wijesinghe, RD Hornung, AL Garcia, NG Hadjiconstantinou, et al.  Three-dimensional hybrid continuum-atomistic simulations for multiscale hydrodynamics.  Transactions of the ASME-I-Journal of Fluids Engineering, 126(5):768–777, 2004.
  * [9] Thomas E Schwartzentruber and Iain D Boyd.  A hybrid particle-continuum method applied to shock waves.  Journal of Computational Physics, 215(2):402–416, 2006.
  * [10] Thomas E Schwartzentruber, Leonardo C Scalabrin, and Iain D Boyd.  A modular particle–continuum numerical method for hypersonic non-equilibrium gas flows.  Journal of Computational Physics, 225(1):1159–1174, 2007.
  * [11] Giacomo Dimarco and Lorenzo Pareschi.  Hybrid multiscale methods ii. kinetic equations.  Multiscale Modeling & Simulation, 6(4):1169–1197, 2008.
  * [12] Jonathan M Burt and Iain D Boyd.  A hybrid particle approach for continuum and rarefied flow simulation.  Journal of Computational Physics, 228(2):460–475, 2009.
  * [13] Pierre Degond, Jian-Guo Liu, and Luc Mieussens.  Macroscopic fluid models with localized kinetic upscaling effects.  Multiscale Modeling & Simulation, 5(3):940–979, 2006.
  * [14] Carlo Cercignani.  The Boltzmann Equation and Its Applications.  Springer, 1988.
  * [15] Kun Xu and Juan-Chen Huang.  A unified gas-kinetic scheme for continuum and rarefied flows.  Journal of Computational Physics, 229(20):7747–7764, 2010.
  * [16] Kun Xu.  Direct modeling for computational fluid dynamics: construction and application of unified gas-kinetic schemes.  World Scientific, 2015.
  * [17] Shi Jin.  Asymptotic preserving (ap) schemes for multiscale kinetic and hyperbolic equations: a review.  Lecture notes for summer school on methods and models of kinetic theory (M&MKT), Porto Ercole (Grosseto, Italy), pages 177–216, 2010.
  * [18] EM Shakhov.  Generalization of the krook kinetic relaxation equation.  Fluid Dynamics, 3(5):95–96, 1968.
  * [19] Tianbai Xiao, Qingdong Cai, and Kun Xu.  A well-balanced unified gas-kinetic scheme for multiscale flow transport under gravitational field.  Journal of Computational Physics, 332:475 – 491, 2017.
  * [20] Clément Mouhot and Lorenzo Pareschi.  Fast algorithms for computing the boltzmann collision operator.  Mathematics of computation, 75(256):1833–1852, 2006.
  * [21] Lei Wu, Craig White, Thomas J Scanlon, Jason M Reese, and Yonghao Zhang.  Deterministic numerical solutions of the boltzmann equation using the fast spectral method.  Journal of Computational Physics, 250:27–52, 2013.
  * [22] Lei Wu, Jason M Reese, and Yonghao Zhang.  Solving the boltzmann equation deterministically by the fast spectral method: application to gas microflows.  Journal of Fluid Mechanics, 746:53–84, 2014.
  * [23] Chang Liu, Kun Xu, Quanhua Sun, and Qingdong Cai.  A unified gas-kinetic scheme for continuum and rarefied flows iv: Full boltzmann and model equations.  Journal of Computational Physics, 314:305–340, 2016.
  * [24] Kun Xu.  A gas-kinetic bgk scheme for the navier–stokes equations and its connection with artificial dissipation and godunov method.  Journal of Computational Physics, 171(1):289–335, 2001.
  * [25] Iain D Boyd, Gang Chen, and Graham V Candler.  Predicting failure of the continuum fluid equations in transitional hypersonic flows.  Physics of Fluids, 7(1):210–219, 1995.
  * [26] Quanhua Sun.  Information preservation methods for modeling micro-scale gas flows.  2003\.
  * [27] Yajun Zhu, Chengwen Zhong, and Kun Xu.  Implicit unified gas-kinetic scheme for steady state solutions in all flow regimes.  Journal of Computational Physics, 315:16–38, 2016.
  * [28] Yajun Zhu, Chengwen Zhong, and Kun Xu.  Unified gas-kinetic scheme with multigrid convergence for rarefied flow study.  Physics of Fluids, 29(9):096102, 2017.

(a) U-velocity (b) Temperature

Figure 22: Density and Temperature contours in the flow passing through
cylinder at $Kn=0.001$ (upper flood: AUGKS, upper lines: UGKS, lower flood:
GKS).

(a) U-velocity (b) Temperature

Figure 23: Density and Temperature contours in the flow passing through
cylinder at $Kn=0.01$ (upper flood: AUGKS, upper lines: UGKS, lower flood:
GKS).

(a) Density (b) U-velocity (c) Temperature

Figure 24: Solutions along the horizontal central line in front of cylinder at
$Kn=0.001$.

(a) Density (b) U-velocity (c) Temperature

Figure 25: Solutions along the horizontal central line behind cylinder at
$Kn=0.001$.

(a) Density (b) U-velocity (c) Temperature

Figure 26: Solutions along the horizontal central line in front of cylinder at
$Kn=0.01$.

(a) Density (b) U-velocity (c) Temperature

Figure 27: Solutions along the horizontal central line behind cylinder at
$Kn=0.01$.

(a) $a_{1}$ (b) $a_{3}$ (c) Collision time (d) Velocity adaptation

Figure 28: Velocity space adaptation in the flow domain at $Kn=0.001$.

(a) $a_{1}$ (b) $a_{3}$ (c) Collision time (d) Velocity adaptation

Figure 29: Velocity space adaptation in the flow domain at $Kn=0.01$.