2023

[1]Peichao Li

1] School of Biomedical Engineering & Imaging Sciences, King’s College London,
London, UK

2] Department of Neurosurgery, King’s College Hospital NHS Foundation Trust,
London, UK

# Spatial gradient consistency for unsupervised learning of hyperspectral
demosaicking: Application to surgical imaging

<EMAIL_ADDRESS>Muhammad Asad  Conor Horgan  Oscar MacCormac  Jonathan
Shapey  Tom Vercauteren  [  [

###### Abstract

Purpose: Hyperspectral imaging has the potential to improve intraoperative
decision making if tissue characterisation is performed in real-time and with
high-resolution. Hyperspectral snapshot mosaic sensors offer a promising
approach due to their fast acquisition speed and compact size. However, a
demosaicking algorithm is required to fully recover the spatial and spectral
information of the snapshot images. Most state-of-the-art demosaicking
algorithms require ground-truth training data with paired snapshot and high-
resolution hyperspectral images, but such imagery pairs with the exact same
scene are physically impossible to acquire in intraoperative settings. In this
work, we present a fully unsupervised hyperspectral image demosaicking
algorithm which only requires exemplar snapshot images for training purposes.

Methods: We regard hyperspectral demosaicking as an ill-posed linear inverse
problem which we solve using a deep neural network. We take advantage of the
spectral correlation occurring in natural scenes to design a novel inter
spectral band regularisation term based on spatial gradient consistency. By
combining our proposed term with standard regularisation techniques and
exploiting a standard data fidelity term, we obtain an unsupervised loss
function for training deep neural networks, which allows us to achieve real-
time hyperspectral image demosaicking.

Results: Quantitative results on hyperspetral image datasets show that our
unsupervised demosaicking approach can achieve similar performance to its
supervised counter-part, and significantly outperform linear demosaicking. A
qualitative user study on real snapshot hyperspectral surgical images confirms
the results from the quantitative analysis.

Conclusion: Our results suggest that the proposed unsupervised algorithm can
achieve promising hyperspectral demosaicking in real-time thus advancing the
suitability of the modality for intraoperative use.

###### keywords:

hyperspectral imaging, demosaicking, unsupervised learning, surgical imaging

##  1 Introduction

Hyperspectral Imaging (HSI) is a technique that captures and processes
spectral data distributed across a large number of wavelengths. It provides a
non-contact, non-ionising and non-invasive solution suitable for many medical
applications (lu2014medical, ; shapey2019intraoperative, ; clancy2020surgical,
). HSI can provide information beyond what human vision can observe, such as
tissue perfusion, oxygen saturation, and other diagnostic measurements
(holmer2018hyperspectral, ). Hence, it facilitates important medical tasks
such as tissue differentiation and characterisation. Depending on the number
of bands, hyperspectral imaging may also be called multispectral imaging, but
in this work we will refer to hyperspectral imaging for consistency.

Snapshot hyperspectral imaging is a promising technique which can capture
hyperspectral images in real-time. Snapshot mosaic cameras are a common type
of snapshot hyperspectral camera which employ multi-spectral filter array
(MSFA) to acquire multi-spectral data in a single exposure. In MSFA cameras
the $n\times n$ sensor arrays are arranged in a repeating pattern similar to
the $2\times 2$ Bayer filter arrays on RGB cameras (Figure 1, left) and are
thus capable of obtaining a maximum of $n^{2}$ bands instantly. However, it
achieves real-time multi-spectral data acquisition at the cost of reducing
both spatial and spectral resolution. Efficient hyperspectral demosaicking
algorithms are thus required to fully restore the spatial and spectral
resolution from the snapshot images. More details on hyperspectral imaging
techniques and snapshot mosaic imaging can be found in (li2022deep, ).

Traditionally, demosaicking algorithms were developed using interpolation-
based methods or statistics-based techniques (yu2006colour, ;
eismann2004application, ), but these methods may still suffer from colour
artifacts and blurriness. Recent deep-learning based algorithms have been
developed for efficient and accurate image super-resolution and demosaicking
tasks. Deep neural networks such as SRCNN (dong2014learning, ), EDSR
(lim2017edsr, ) and RNAN (zhang2019residual, ) have demonstrated their
performances on RGB image super-resolution tasks, and thus similar methods
have been extended to process hyperspectral images (mei2017hyperspectral, ;
dijkstra2019hyperspectral, ). (arad2022ntire, ) introduced several state-of-
the-art learning-based hyperspectral demosaicking algorithms of natural scenes
in NTIRE 2022 Spectral Demosaicking Challenge. The leading contestants include
Enhanced HAN (niu2020single, ), NLRAN (arad2022ntire, ) and Res2-Unet based
methods (song2022hyperspectral, ). Our previous work (li2022deep, ) also
demonstrated the use of a synthetic surgical HSI dataset and deep-learning
models for developing hyperspectral demosaicking algorithms suitable for
intraoperative surgical guidance tasks.

However, most deep-learning based demosaicking algorithms rely on a large
number of high-resolution HSI data as the ground truth for model training.
Publicly available medical hyperspectral datasets such as HELICoiD
(fabelo2019helicoid, ) and ODSI (hyttinen2020oral, ) involve large line-scan
or spectral-scan HSI systems to obtain high-resolution hyperspectral data, and
the acquisition speed is slow. Consequently, these imaging systems are not
ideal for intraoperative use. Fortunately, (ebner2021intraoperative, )
demonstrated that the acquisition of intraoperative snapshot mosaic images is
less challenging as its compact imaging system can be seamlessly integrated
into a standard surgical workflow.

This paper presents an unsupervised-learning-based HSI demosaicking algorithm
which uses only snapshot mosaic images and does not require corresponding
high-resolution images for training. A demosaicking loss function is proposed
based on a novel spatial gradient consistency regularisation technique
combined with traditional regularisation methods including Tikhonov
regularisation and total variation. The proposed algorithm has been tested
with 3 different deep neural networks on 3 different datasets. Quantitative
measures have been performed to compare the unsupervised algorithm against
linear demosaicking and supervised training, and a qualitative user study was
conducted to validate the proposed algorithm on a medical HSI dataset.

Figure 1:  (Left) $4\times 4$ MSFA for an IMEC snapshot camera, the colour of
each pixel correlates to the perceived colour of a human observer. (Middle)
Spectral responses of all 16 sensors on the MSFA of an IMEC snapshot camera.
(Right) Wasserstein metric heatmap measuring distances between different
spectral responses of the sensors on an IMEC snapshot camera.

##  2 Materials and methods

###  2.1 Demosaicking as an ill-posed linear inverse problem

Problem formulation. Hyperspectral image demosaicking involves recovering the
fully sampled hyperspectral image $I\in\mathbb{R}^{X\times Y\times C}$ from a
snapshot image $I^{s}\in\mathbb{R}^{X\times Y}$, where $X$ and $Y$ are the
spatial dimensions and $C$ is the number of spectral bands. The relationship
between $I$ and $I^{s}$ can be expressed through a linear degradation operator
$\mathcal{D}$:

$I^{s}=\mathcal{D}(I)$ (1)

For a typical MSFA arrangement as shown in Figure 1 (left), $\mathcal{D}$ can
be simply expressed as a selection matrix containing only 0 and 1, thereby
mapping the pixel values of $I^{s}$ from $I$. In other words, for each spatial
location $(x,y)$, there is a single corresponding spectral band $c_{x,y}$ such
that $I^{s}(x,y)=I(x,y,c_{x,y})$ The inverse problem corresponding to (1) is
ill-posed because of the highly ill-conditioned selection operator
$\mathcal{D}$. Therefore, appropriate regularisation is required. A classical
inverse problem approach would aim at solving for

$\hat{I}=\arg\min_{I}\big{[}\mathcal{L}(I^{s},\mathcal{D}(I))+\lambda\mathcal{R}(I)\big{]}$
(2)

where $\mathcal{L}(I^{s},\mathcal{D}(I))$ is the data fidelity term that
measures the differences between the the known snapshot image $I^{s}$ and the
subsampling of the unknown fully-sampled hyperspectral image $I$.
$\mathcal{R}$ represents the regularisation terms. $\lambda$ is the
regularisation factor that determines the trade-off between the data fidelity
and regularisation.

Translating this into an unsupervised machine learning setting, we now seek to
optimise for the parameters $\theta$ of a deep neural network $f_{\theta}$
mapping a snapshot mosaic input $I^{s}$ to a fully-sampled hyperspectral image
$f_{\theta}(I^{s})$:

$\hat{\theta}=\arg\min_{\theta}\mathbb{E}_{I^{s}}\big{[}\mathcal{L}(I^{s},\mathcal{D}(f_{\theta}(I^{s})))+\lambda\mathcal{R}(f_{\theta}(I^{s}))\big{]}$
(3)

where the expectation $\mathbb{E}_{I^{s}}$ is to be considered as being taken
over an empirical distribution defined by a training set of snapshot mosaic
images (with no need for ground truth).

Spatial gradient consistency regularisation. Regularisation terms in (3) aim
at incorporating prior information about the problem being solved. In our
case, all spectral bands are imaging the same physical scene. We also observe
that the spectrum of natural objects and biological tissues present with
specific characteristics such as continuity and smoothness. Additionally, the
response functions corresponding to the different spectral bands as shown in
Figure 1 (middle) shares significant spectral overlap. It is thus expected
that our spectral bands will exhibit substantial correlation. Inter-spectral
band correlation was notably demonstrated empirically for RGB images in
(gunturk2002color, ). However, while correlation is expected, assuming a
simple linear relationship would make for too crude an approximation.

Here, inspired by image similarity metrics that exploit image gradients for
multimodal image registration where non-trivial correlation across the imaging
modalities is expected (haber2006intensity, ), we propose to promote
correlation between the spatial gradients of the individual spectral bands in
our reconstructions. Let $c_{1}$ and $c_{2}$ be the indices of two spectral
bands of interest, with $I^{c}=I(\cdot,\cdot,c)$, and $c\in(c_{1},c_{2})$ the
corresponding spectral band images. For simplicity, we make use of forward
differences to compute spatial gradients:
$\nabla_{x}I^{c}(x,y)=I^{c}(x+1,y)-I^{c}(x,y)$ and
$\nabla_{y}I^{c}(x,y)=I^{c}(x,y+1)-I^{c}(x,y)$. We propose to consider the
correlation coefficient between the spatial gradients as a regularisation:

$\mathcal{R}^{c_{1},c_{2}}_{\operatorname{\rho}}(I)=-\operatorname{\rho}(\nabla_{x}I^{c_{1}},\nabla_{x}I^{c_{2}})-\operatorname{\rho}(\nabla_{y}I^{c_{1}},\nabla_{y}I^{c_{2}})$
(4)

Given $C$ spectral bands, $C^{2}$ pairwise comparisons are possible. However,
the strength of the correlation is not expected to be the same for all pairs
of bands. Indeed, two bands with close spectral peaks should lead to higher
correlation than two bands with further peaks. Given the complex structure of
the spectral response functions shown in Figure 1 (middle), we propose to
weight the contribution of each pair of spectral band according to the
Wasserstein distance $W_{c_{1},c_{2}}$ between the spectral response functions
of the two bands:

$\mathcal{R}_{\operatorname{\rho}}(I)=\sum_{c_{1}\neq
c_{2}}e^{-\frac{W_{c_{1},c_{2}}}{\tau}}\leavevmode\nobreak\
\mathcal{R}^{c_{1},c_{2}}_{\operatorname{\rho}}(I)$ (5)

where the negative exponential mapping with temperature scaling $\tau$ allows
to control the relative importance of each pair. The exponential Wasserstein
distance gives an indication of how closely the spectral responses of the two
bands might be correlated, as shown in the heatmap in Figure 1 (right), where
lighter colour means the two spectral bands are closer. By strengthening the
correlation between the spatial gradient maps of different spectral bands we
expect to enhance the sharp edges and contours.

Other regularisation terms. Tikhonov regularisation is a common method for
ill-conditioned problems. It can be characterised as:

$\mathcal{R}_{\textrm{Tik}}(I)=\|\bm{\Gamma}\cdot I\|_{2}^{2}$ (6)

Here, we choose to use the Laplacian matrix as the Tikhonov matrix
$\bm{\Gamma}$ to deal with potential high-frequency artifacts introduced
during the super-resolution process. While Tikhonov regularisation can
effectively eliminate undesirable outliers and led to smooth images, it also
has the potential risk of applying too much smoothness and erasing all sharp
edges and contours, which is harmful for recovering details in the images.

Total variation is another term which is able to preserve edges while
regularising solutions of the inverse problem:

$\mathcal{R}_{\textrm{TV}}(I)=\|\nabla_{x}I\|_{1}+\|\nabla_{y}I\|_{1}$ (7)

By combining our proposed spatial gradient consistency term with Tikhonov and
total variation regularisation, we obtain the regularisation term
$\mathcal{R}$ in (2) using $\lambda_{\textrm{Tik}}$, $\lambda_{\textrm{TV}}$
and $\lambda_{\operatorname{\rho}}$ as weighting factors for individual terms:

$\mathcal{R}(I)=\lambda_{\textrm{Tik}}\mathcal{R}_{\textrm{Tik}}(I)+\lambda_{\textrm{TV}}\mathcal{R}_{\textrm{TV}}(I)+\lambda_{\operatorname{\rho}}\mathcal{R}_{\operatorname{\rho}}(I)$
(8)

###  2.2 Image demosaicking pipeline

Figure 2: The pipeline of the proposed unsupervised demosaicking algorithm.

Figure 2 depicts the general pipeline of our proposed algorithm using deep
neural networks for hyperspectral image demosaicking problems. It starts from
the input snapshot mosaic images where bilinear interpolation based
demosaicking can be applied to recover the spatial and spectral dimension of
the images. The linearly interpolated images serve as the input of the network
to generate refined demosaicking results. Most deep neural networks for image
super-resolution or demosaicking can be integrated into this pipeline, such as
U-Net (ronneberger2015unet, ), EDSR (lim2017edsr, ) and Res2-Unet
(song2022hyperspectral, ).

Aside from the network, given that the measured pixels in the original
snapshot $I^{s}$ should be equal to the corresponding pixels in the
demosaicked hypercube $I$, we propose to include an overriding operator which
applies the pixel values from $I^{s}$ to their corresponding position in $I$.
This forces the data fidelity term $\mathcal{L}$ in (2) to be always 0
irrespective of the metric we choose. Based on the output images from the
network with the overridden snapshot pixels, the Tikhonov regularisation,
total variation and the spatial gradient consistency regularisation terms are
calculated and minimised using gradient descent, and the parameters in the
networks are updated.

###  2.3 Source datasets

To experiment the proposed demosaicking algorithm, three hyperspectral imaging
datasets are used in this work, which will be presented in this section.

HELICoiD. (fabelo2019helicoid, ) presented a publicly available in-vivo
hyperspectral human brain image dataset within the European project HELICoiD
(HypErspectraL Imaging Cancer Detection). The hyperspectral images in this
dataset were acquired using a line-scan hyperspectral camera system capable of
capturing high spectral-resolution hypercubes during neurosurgical operations.
The dataset contains 36 images in the Visual and Near Infrared (VNIR) range
from 400nm to 1000nm. We applied the same method described in (li2022deep, )
to perform white balancing, and then simulated snapshot mosaic images and
their corresponding high-resolution demosaicked hypercubes using spectral
response functions of a real hyperspectral snapshot camera.

ARAD_1K. With the NTIRE 2022 Spectral Demosaicking Challenge, (arad2022ntire,
) provided 1000 hyperspectral images of natural scenes with 16 spectral bands
ranging from 400nm to 1000nm. The snapshot images were simulated following a
$4\times 4$ MSFA pattern. There were 950 hyperspectral images for training,
where the simulated snapshot images and their corresponding ground truth
images were both provided. The other 50 images were for testing, but the
ground truth was not publicly available, so we separated 50 images out from
the 950 training set for testing.

NeuroHSI. NeuroHSI is an actively running, NIHR funded, single centre
prospective observational study assessing the intra-operative capabilities of
a $4\times 4$, 16 band visible range snapshot mosaic camera (IMEC
CMV2K-SSM4X4-VIS) to differentiate between pathological tissue and healthy
brain tissue, as well as to evaluate custom made algorithms capable of
correlating information from specific bands to tissue oxygenation
measurements. Phase one of this study has now been completed and video
hyperspectral data from two brain metastases, two gliomas (WHO grades 2-4),
one meningiomas, one vestibular schwannoma, one cerebral aneurysm and one
cerebral arteriovenous malformation has been collected. 150 snapshot images
with minor motion blur or out-of-focus blur were manually selected from the
video data of the 8 patients, where 90 images from 4 patients are reserved for
training, 30 images from 2 patients reserved for validation and 30 from the
remaining 2 patients for testing.

###  2.4 Implementation details

Our proposed algorithm was implemented with PyTorch and tested on all three
datasets described in Section 2.3. For the HELICoiD dataset, synthetic
snapshot images and their corresponding high-resolution hypercubes were
simulated using sensor information from the snapshot camera IMEC
CMV2K-SSM4X4-VIS. The dataset was divided into 3 groups: 24 images acquired
from 15 different patients as the training set, 6 images from 4 patients as
the validation set, and the remaining 6 images from 3 patients as the test
set. For the ARAD_1K dataset, the original raw snapshot data were simulated
with an unknown exposure setting. Recovering such an unknown exposure is not
the primary focus for our experiment. Therefore, new snapshot images were
simulated using the ground truth hypercubes and the MSFA simulation algorithm
provided by the organiser. The dataset was also divided into 3 groups: 720
images for training, 180 for validation and 50 for testing.

As both the HELICoiD and ARAD_1K datasets have high-resolution hypercubes as
ground truths, the U-Net, EDSR and Res2-Unet models were trained in both a
supervised and an unsupervised manner. For supervised training, the models
were all trained using the Mean Relative Absolute Error (MRAE) Loss as
described in (song2022hyperspectral, ). For unsupervised training, the
regularisation terms described in (8) were used as the loss function, and the
models were trained with only the simulated snapshot images as inputs. The
regularisation factors in (8) were set to $\lambda_{\textrm{Tik}}=1$,
$\lambda_{\textrm{TV}}=10^{-3}$ and $\lambda_{\operatorname{\rho}}=1$
respectively, and the temperature scaling $\tau$ in (5) was set to 0.1.
Details on the parameter selection and the ablation study can be found in the
supplementary material. Random flipping and rotation were not performed
because they can disrupt the MSFA pattern of the snapshot images. Therefore,
random divisible spatial cropping were performed where the position and size
of the crop were all divisible by the size of the mosaic. The network models
were trained using the Adam optimiser with $\beta_{1}=0.5$ and
$\beta_{2}=0.99$ and a batch size of 4. The initial learning rate was set to
$1\times 10^{-4}$. Results were quantitatively evaluated based on 3 metrics,
including Structural Similarity (SSIM), Peak Signal-to-Noise Ratio (PSNR) and
Spectral Angle Mapper (SAM) (kruse1993spectral, ).

The 150 image frames selected from the NeuroHSI video dataset were all
acquired from an IMEC CMV2K-SSM4X4-VIS camera, and there are no ground truth
high-resolution hypercubes, so the experiment only involves unsupervised
training. 90 snapshot image frames from 4 patients were used for training, and
30 images from 2 patients for both validation and testing. Res2-Unet was
adopted for the proposed algorithm, and the parameters used for training on
NeuroHSI dataset remains the same as the HELICoiD and ARAD_1K dataset. The
results were evaluated qualitatively by a user study which will be described
in Section 3.2.

##  3 Results

Table 1: Comparison of demosaicking accuracy between linear demosaicking and
different networks with supervised and unsupervised training setup on HELICoiD
and ARAD_1K datasets.

Dataset | Method | Network | SSIM | PSNR | SAM
---|---|---|---|---|---
HELICoiD | Linear | - | $0.969\pm 0.011$ | $32.5\pm 2.48$ | $0.123\pm 0.051$
Supervised | UNet | $0.991\pm 0.005$ | $38.7\pm 3.06$ | $0.057\pm 0.025$
EDSR | $0.994\pm 0.004$ | $42.3\pm 3.70$ | $0.035\pm 0.013$
Res2-Unet | $\mathbf{0.995\pm 0.004}$ | $\mathbf{42.9\pm 3.99}$ | $\mathbf{0.032\pm 0.012}$
Unsupervised | UNet | $0.991\pm 0.004$ | $38.6\pm 3.43$ | $0.057\pm 0.024$
EDSR | $\mathbf{0.993\pm 0.004}$ | $\mathbf{39.9\pm 3.70}$ | $\mathbf{0.046\pm 0.019}$
Res2-Unet | $0.992\pm 0.005$ | $39.7\pm 3.66$ | $0.051\pm 0.019$
ARAD_1K | Linear | - | $0.956\pm 0.039$ | $33.4\pm 4.16$ | $0.055\pm 0.022$
Supervised | UNet | $0.978\pm 0.022$ | $38.4\pm 3.91$ | $0.041\pm 0.016$
EDSR | $0.993\pm 0.007$ | $43.1\pm 3.74$ | $0.025\pm 0.009$
Res2-Unet | $\mathbf{0.998\pm 0.002}$ | $\mathbf{48.7\pm 3.88}$ | $\mathbf{0.015\pm 0.006}$
Unsupervised | UNet | $0.977\pm 0.020$ | $38.1\pm 3.62$ | $0.041\pm 0.015$
EDSR | $\mathbf{0.989\pm 0.013}$ | $41.0\pm 4.07$ | $0.031\pm 0.013$
Res2-Unet | $0.988\pm 0.017$ | $\mathbf{41.3\pm 4.42}$ | $\mathbf{0.031\pm 0.012}$
|  |  |  |  |

###  3.1 Quantitative evaluation

The quantitative results of the demosaicked hypercubes on both HELICoiD and
ARAD_1K datasets are shown in Table 1. Paired T-test was performed to compare
against the performance of two demosaicking methods. For both datasets, the
supervised training of Res2-Unet achieved the highest demosaicking accuracy.
The supervised EDSR results did not show statistical differences compared to
Res2-Unet at a significant level of 0.05 on the HELICoiD dataset, with
p-values of 0.35, 0.34 and 0.30 for SSIM, PSNR and SAM respectively. However,
on the ARAD_1K dataset the p-values of $<10^{-5}$ for all 3 metrics indicates
that Res2-Unet outperforms EDSR significantly.

The demosaicking results of the proposed unsupervised method on Res2-Unet are
significantly lower than the supervised method with p-values of 0.040, 0.016,
0.007 on the 3 metrics on HELICoiD dataset, and p-values of close to 0 on
ARAD_1K dataset, showing that our proposed method cannot match state-of-the-
art supervised demosaicking methods when ground truths are provided. However,
when comparing supervised and unsupervised EDSR results, the p-values of 0.17,
0.06 and 0.07 on the HELICoiD dataset indicates that our proposed method can
still reach similar performance as a supervised method. On the ARAD_1K
dataset, although the unsupervised EDSR performs significantly lower than
supervised EDSR with p-values of 0.02, 0.0001 and 0.0005, it still outperforms
the supervised U-Net significantly with p-values of $<10^{-5}$ for all 3
metrics. In both datasets, all supervised and unsupervised results
significantly outperform linear demosaicking with p-values close to 0.

The speed of our proposed demosaicking algorithm depends on the choice of
network. For a single image of size $512\times 480$ from the ARAD_1K dataset,
the inference times for UNet, EDSR and Res2-Unet are around 0.009s, 0.006s and
0.010s respectively with NVIDIA RTX 3080 Ti. This demonstrates that when
combining a suitable neural network and computing hardware, our proposed
algorithm can achieve high quality hyperspectral demosaicking in real-time.

Figure 3: Comparison between different demosaicking methods on an example
NeuroHSI test image. The reconstructed sRGB images are converted from the
demosaicked hyperspectral data following the method described in (li2022deep,
).

###  3.2 Qualitative evaluation and user study

As there is no ground truth data for the NeuroHSI dataset, a qualitative user
study was conducted to evaluate the demosaicked results of the NeuroHSI
dataset. The user study was conducted using forced-choice pairwise comparison
(mantiuk2012image, ). Figure 3 illustrates the pseudo-sRGB reconstructions of
an example NeuroHSI patient image tested using three methods: linear
demosaicking (L), supervised Res2-Unet model trained from HELICoiD dataset
(SL) and the unsupervised Res2-Unet model trained from NeuroHSI training set
(UL). 30 test images were included in the user study, each tested with the
three methods (L, SL, UL). There are thus 90 questions in total, each
containing two images of the same scene with 2 different demosaicking methods.
These questions were divided into 3 separate surveys, each containing 30
questions. Participants were randomly assigned to answer one of 3 surveys and
asked to choose the image with better quality for each question (pair of
images) without any knowledge of which demosaicking method was used. The
participants of this survey were all neurosurgical experts with 2 to 15 years
of experience. We received 12 responses in total that are summarised in Table
2. We applied the Bradley-Terry model (bradley1952rank, ) to rank the
demosaicking methods, which gives the estimated preference scale of
$\pi=(0.050,0.445,0.505)$ for L, SL and UL respectively. This indicates that
the experts considered the images recovered from our proposed demosaicking
method to have similar quality as the images from a supervised model, with the
baseline linear demosaicking the least favourable method. More details can be
found in the supplementary material.

Table 2: Number of votes received for each demosaicking method in all pairwise comparisons in the image quality assessment survey. Linear | Supervised | Linear | Unsupervised | Supervised | Unsupervised
---|---|---|---|---|---
13 | 107 | 10 | 110 | 57 | 63
|  |  |  |  |

##  4 Conclusion

In this work, we have presented a novel unsupervised approach for medical
hyperspectral image demosaicking. The proposed algorithm does not rely on
high-resolution medical hyperspectral data which are hard to acquire in a
surgical environment, but instead only snapshot mosaic images are required,
which are much easier to capture. The combination of Tikhonov regularisation,
total variation and spectral correlation regularisation has been adopted for
unsupervised network training, and the results were tested both quantitatively
and qualitatively, showing convincing results over basic linear demosaicking,
and comparable results against supervised demosaicking methods, thus proving
its capability for real-time intraoperative surgical application.

Supplementary information For more information regarding details of the
ablation study, quantitative metrics, qualitative results and the user study,
please refer to the supplementary document alongside this article.

## Declarations

  * •

Funding: This study/project is funded by the NIHR [NIHR202114]. The views
expressed are those of the author(s) and not necessarily those of the NIHR or
the Department of Health and Social Care. This work was supported by core
funding from the Wellcome/EPSRC [WT203148/Z/16/Z; NS/A000049/1]. This project
has received funding from the European Union’s Horizon 2020 research and
innovation programme under grant agreement No 101016985 (FAROS project). TV is
supported by a Medtronic / RAEng Research Chair [RCSRF1819\7\34]. PL is funded
by China Scholarship Council. CH is supported by an InnovateUK Secondment
Scholars Grant (Project Number 75124). For the purpose of open access, the
authors have applied a CC BY public copyright license to any Author Accepted
Manuscript version arising from this submission.

  * •

Conflict of interest: TV and JS are co-founders and shareholders of
Hypervision Surgical.

  * •

Ethics approval: All procedures within this study involving human subjects
were in accordance with both the institutional and regional ethical committee
(REC reference 22/LO/0046, IRAS 284230) and with the 1964 Helsinki declaration
and its later amendments.

  * •

Informed consent: Informed consent was obtained from all individual
participants involved in the study.

  * •

Consent for publication: The authors affirm that human research participants
provided informed consent for publication of the images in Figure 3.

## References

  *   * (1) Lu, G., Fei, B.: Medical hyperspectral imaging: a review. Journal of Biomedical Optics 19(1), 010901 (2014)
  * (2) Shapey, J., Xie, Y., Nabavi, E., Bradford, R., Saeed, S.R., Ourselin, S., Vercauteren, T.: Intraoperative multispectral and hyperspectral label-free imaging: A systematic review of in vivo clinical studies. Journal of Biophotonics 12(9), 201800455 (2019)
  * (3) Clancy, N.T., Jones, G., Maier-Hein, L., Elson, D.S., Stoyanov, D.: Surgical spectral imaging. Medical image analysis 63, 101699 (2020)
  * (4) Holmer, A., Marotz, J., Wahl, P., Dau, M., Kämmerer, P.W.: Hyperspectral imaging in perfusion and wound diagnostics – methods and algorithms for the determination of tissue parameters. Biomedical Engineering / Biomedizinische Technik 63(5), 547–556 (2018). https://doi.org/10.1515/bmt-2017-0155
  * (5) Li, P., Ebner, M., Noonan, P., Horgan, C., Bahl, A., Ourselin, S., Shapey, J., Vercauteren, T.: Deep learning approach for hyperspectral image demosaicking, spectral correction and high-resolution RGB reconstruction. Computer Methods in Biomechanics and Biomedical Engineering: Imaging & Visualization 10(4), 409–417 (2022). https://doi.org/10.1080/21681163.2021.1997646
  * (6) Yu, W.: Colour demosaicking method using adaptive cubic convolution interpolation with sequential averaging. Vision, Image and Signal Processing 153, 666–676 (2006). https://doi.org/10.1049/ip-vis:20050281
  * (7) Eismann, M.T., Hardie, R.C.: Application of the stochastic mixing model to hyperspectral resolution enhancement. IEEE Transactions on Geoscience and Remote Sensing 42(9), 1924–1933 (2004). https://doi.org/10.1109/TGRS.2004.830644
  * (8) Dong, C., Loy, C.C., He, K., Tang, X.: Learning a deep convolutional network for image super-resolution. In: Fleet, D., Pajdla, T., Schiele, B., Tuytelaars, T. (eds.) ECCV 2014, pp. 184–199. Springer, Cham (2014)
  * (9) Lim, B., Son, S., Kim, H., Nah, S., Lee, K.M.: Enhanced deep residual networks for single image super-resolution. In: CVPR Workshops 2017, pp. 1132–1140 (2017). https://doi.org/10.1109/CVPRW.2017.151
  * (10) Zhang, Y., Li, K., Li, K., Zhong, B., Fu, Y.: Residual non-local attention networks for image restoration. In: ICLR 2019. OpenReview, ??? (2019)
  * (11) Mei, S., Yuan, X., Ji, J., Zhang, Y., Wan, S., Du, Q.: Hyperspectral image spatial super-resolution via 3d full convolutional neural network. Remote Sensing 9(11) (2017). https://doi.org/10.3390/rs9111139
  * (12) Dijkstra, K., van de Loosdrecht, J., Schomaker, L., Wiering, M.A.: Hyperspectral demosaicking and crosstalk correction using deep learning. Machine Vision and Applications 30(1), 1–21 (2019)
  * (13) Arad, B., Timofte, R., Yahel, R., Morag, N., Bernat, A., Wu, Y., Wu, X., Fan, Z., Xia, C., Zhang, F., Liu, S., Li, Y., Feng, C., Lei, L., Zhang, M., Feng, K., Zhang, X., Yao, J., Zhao, Y., Ma, S., He, F., Dong, Y., Yu, S., Qiu, D., Liu, J., Bi, M., Song, B., Sun, W., Zheng, J., Zhao, B., Cao, Y., Yang, J., Cao, Y., Kong, X., Yu, J., Xue, Y., Xie, Z.: NTIRE 2022 spectral demosaicing challenge and data set. In: 2022 IEEE/CVF CVPR Workshops, pp. 881–895 (2022). https://doi.org/10.1109/CVPRW56347.2022.00103
  * (14) Niu, B., Wen, W., Ren, W., Zhang, X., Yang, L., Wang, S., Zhang, K., Cao, X., Shen, H.: Single image super-resolution via a holistic attention network. In: Vedaldi, A., Bischof, H., Brox, T., Frahm, J.-M. (eds.) ECCV 2020, pp. 191–207. Springer, Cham (2020)
  * (15) Song, B., Ma, S., He, F., Sun, W.: Hyperspectral reconstruction from RGB images based on Res2-Unet deep learning network. Optics and Precision Engineering 30(13), 1606 (2022)
  * (16) Fabelo, H., Ortega, S., Szolna, A., Bulters, D., Piñeiro, J.F., Kabwama, S., J-O’Shanahan, A., Bulstrode, H., Bisshopp, S., Kiran, B.R., Ravi, D., Lazcano, R., Madroñal, D., Sosa, C., Espino, C., Marquez, M., De La Luz Plaza, M., Camacho, R., Carrera, D., Hernández, M., Callicó, G.M., Morera Molina, J., Stanciulescu, B., Yang, G.-Z., Salvador, R., Juárez, E., Sanz, C., Sarmiento, R.: In-vivo hyperspectral human brain image database for brain cancer detection. IEEE Access 7, 39098–39116 (2019). https://doi.org/10.1109/ACCESS.2019.2904788
  * (17) Hyttinen, J., Fält, P., Jäsberg, H., Kullaa, A., Hauta-Kasari, M.: Oral and dental spectral image database—odsi-db. Applied Sciences 10(20) (2020). https://doi.org/10.3390/app10207246
  * (18) Ebner, M., Nabavi, E., Shapey, J., Xie, Y., Liebmann, F., Spirig, J.M., Hoch, A., Farshad, M., Saeed, S.R., Bradford, R., Yardley, I., Ourselin, S., Edwards, A.D., Führnstahl, P., Vercauteren, T.: Intraoperative hyperspectral label-free imaging: from system design to first-in-patient translation. Journal of Physics D: Applied Physics 54(29), 294003 (2021). https://doi.org/10.1088/1361-6463/abfbf6
  * (19) Gunturk, B.K., Altunbasak, Y., Mersereau, R.M.: Color plane interpolation using alternating projections. IEEE Transactions on Image Processing 11(9), 997–1013 (2002). https://doi.org/10.1109/TIP.2002.801121
  * (20) Haber, E., Modersitzki, J.: Intensity gradient based registration and fusion of multi-modal images. In: Larsen, R., Nielsen, M., Sporring, J. (eds.) Medical Image Computing and Computer-Assisted Intervention – MICCAI 2006, pp. 726–733. Springer, Berlin, Heidelberg (2006)
  * (21) Ronneberger, O., Fischer, P., Brox, T.: U-Net: Convolutional networks for biomedical image segmentation. In: MICCAI 2015, pp. 234–241 (2015). Springer
  * (22) Kruse, F.A., Lefkoff, A.B., Boardman, J.W., Heidebrecht, K.B., Shapiro, A., Barloon, P.J., Goetz, A.F.H.: The spectral image processing system (sips) interactive visualization and analysis of imaging spectrometer data. Remote Sensing of Environment 44, 145–163 (1993)
  * (23) Mantiuk, R.K., Tomaszewska, A., Mantiuk, R.: Comparison of four subjective methods for image quality assessment. Computer Graphics Forum 31(8), 2478–2491 (2012). https://doi.org/%****␣combined.bbl␣Line␣500␣****10.1111/j.1467-8659.2012.03188.x
  * (24) Bradley, R.A., Terry, M.E.: Rank analysis of incomplete block designs: I. the method of paired comparisons. Biometrika 39(3/4), 324–345 (1952)

See pages - of supplementary.pdf


# Elastic stability of Ga2O3: Addressing the $\beta$ to $\alpha$ phase
transition from first principles

Konstantin Lion  Pasquale Pavone  Claudia Draxl  Institut für Physik and IRIS
Adlershof, Humboldt-Universität zu Berlin, Berlin, Germany

###### Abstract

Elastic and structural properties of $\beta$-Ga2O3 and $\alpha$-Ga2O3 are
investigated from first principles. The full elastic tensors and elastic
moduli of both phases at $0\text{\,}\mathrm{K}$ are computed in the framework
of semi-local density-functional theory. We determine mechanical instabilities
of $\beta$-Ga2O3 by evaluating the full stiffness tensor under load for a
range of hydrostatic pressure values. While a phase transition from the
$\beta$ to $\alpha$ phase is found to be energetically favored at
$2.6\text{\,}\mathrm{GPa}$, we show that the $\beta$ phase is only
mechanically unstable for much higher pressures
($>$30\text{\,}\mathrm{GPa}$$), which agrees well with experimental results.
Our employed approach is based on the Born stability criterion, is independent
of crystal symmetry, and thus can be readily applied to different materials.

##  I Introduction

The wide-gap transparent conducting oxide gallium oxide, Ga2O3, has gained a
lot of interest in recent years as a potential candidate for a number of
applications. Its tunable electrical and optical properties make it a
promising material for gas sensors [1, 2, 3, 4, 5], field-effect transistors
[6], and photodetectors [7, 8, 9, 10]. The material exhibits polymorphism,
i.e., depending on the experimental conditions, it can adopt one out of at
least five different known structures ($\alpha$, $\beta$, $\gamma$, $\delta$,
and $\epsilon$) [11]. The thermodynamically stable phase at ambient conditions
is $\beta$-Ga2O3. It crystallizes in a base-centered monoclinic structure
(space-group C$2/m$) and consists of both tetrahedrally and octahedrally
coordinated gallium atoms. The metastable high-density $\alpha$ phase exhibits
a rhombohedral corundum structure (space-group $R\bar{3}c$) and is solely made
up of octahedrally coordinated gallium. The unit cells of $\alpha$\- and
$\beta$-Ga2O3 are illustrated in Fig. 1. Remeika and Marezio first reported a
phase transition of the $\beta$ phase to the $\alpha$ phase at
$4.4\text{\,}\mathrm{GPa}$ and $1000\text{\,}\mathrm{K}$ [12]. The
transformation was found to be irreversible after quenching the sample to room
temperature. Since then, the phase transition has been the subject of several
studies at both high and room temperatures. Studying nanoparticles of
$\beta$-Ga2O3 embedded in an amorphous silica-based host matrix by
synchrotron-radiation-based x-ray diffraction [13], Lipinska-Kalita and
coworkers reported a phase transition that sets in at about
$6\text{\,}\mathrm{GPa}$ and is not completed at $15\text{\,}\mathrm{GPa}$. It
is, however, not clear whether the phase transition is affected by the host
matrix or intrinsic to the nanoparticles [13]. In a follow-up study conducted
on bulk $\beta$-Ga2O3 at pressures up to $70\text{\,}\mathrm{GPa}$, they
reported an onset of $7.9\text{\,}\mathrm{GPa}$ ($3\text{\,}\mathrm{GPa}$) and
the completion of the transition at $40\text{\,}\mathrm{GPa}$
($30\text{\,}\mathrm{GPa}$) with (without) nitrogen as a pressure-transmitting
medium [14]. Wang and coworkers subsequently subjected freestanding
$\beta$-Ga2O3 nanocrystals to pressures up to $64.9\text{\,}\mathrm{GPa}$ at
room temperature and reported a transition onset between $13.6$ and
$16.4\text{\,}\mathrm{GPa}$ and a completed transition at
$39.2\text{\,}\mathrm{GPa}$ [15]. A phase transition was also observed in
$\beta$-Ga2O3 microparticles, taking place between 20 and
$39\text{\,}\mathrm{GPa}$, where only a highly disordered structure comparable
to $\alpha$-Ga2O3 remained [16]. Also in recently reported shock-recovery
experiments, a phase transition occurred (between $11$ and
$16\text{\,}\mathrm{GPa}$) [17]. Similar transition pressures have been
published in computational studies. The reported values range from $2$ to
$17\text{\,}\mathrm{GPa}$ [18, 19, 20, 21].

Figure 1: Top: Conventional (left) and primitive (right) unit cell of
$\beta$-Ga2O3. Bottom: Same for $\alpha$-Ga2O3. Gallium atoms are in green,
oxygen atoms in red.

From the review above, it is apparent that while there is an extensive
discussion in the literature about the phase transition between the two
polymorphs, it is far from being settled at which pressure the transition
takes place. In particular, the interplay between the transition pressure
$p_{\text{t}}$, that is obtained by thermodynamical considerations, and the
critical pressure $p_{\text{c}}$, that is ruled by a mechanical instability,
has not been thoroughly discussed. This is the aim of our work, where we
illustrate that only considering both allows for understanding the diversity
of experimental results. We investigate the $\beta$ to $\alpha$ phase
transition by calculating the elastic properties of both $\alpha$\- and
$\beta$-Ga2O3, and examining, in particular, their variation under hydrostatic
pressure. The critical pressure is obtained by applying the generalized Born
stability criterion to the elastic constants under load.

##  II Theoretical background

Application of strain leads to a deformation of the crystal due to the
resulting stress. In the linear elastic regime, strain and stress can be
represented by symmetric second-rank tensors and are related by Hooke’s law
[22, 23, 24]

$\sigma_{ij}=C_{ijkl}\,\epsilon_{kl}\,,$ (1)

where $C_{ijkl}$ represents the stiffness tensor, $\epsilon_{kl}$ the physical
strain tensor, and $\sigma_{ij}$ the stress tensor. For conciseness, we employ
Einstein’s notation for summations throughout this work. When dealing with
symmetric tensors, it is convenient to use the Voigt notation, where a pair of
Cartesian indices $ij$ can be written as a single index $\gamma$, according to

$\displaystyle\begin{array}[]{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}ij&11&22&33&23&13&12\\\
\gamma&1&2&3&4&5&6\end{array}.$

In the following, tensors expressed in the Voigt notation are represented by
six-dimensional vectors using a bold font. In order to calculate the stiffness
tensor, one can expand the total crystal energy per unit cell,
$E(\bm{\epsilon})$, in terms of the strain [25]

$\frac{E(\bm{\epsilon})}{V_{0}}=\frac{E_{0}}{V_{0}}+\underbrace{\bm{\sigma}_{0}^{\phantom{I}}\cdot\bm{\epsilon}}_{\text{{\hbox{=}}
{\hbox{0}}}}+\frac{1}{2!}\,\bm{\epsilon}^{\intercal}\cdot\bm{C}\cdot\bm{\epsilon}+\dotsc,$
(2)

where $E_{0}$ and $V_{0}$ refer to the equilibrium energy and volume per unit
cell and $\bm{\sigma}_{0}^{\phantom{I}}$ to the equilibrium stress tensor. The
stiffness tensor is defined as

$C_{\gamma\lambda}=\frac{1}{V_{0}}\,\frac{\partial^{2}\,E(\bm{\epsilon})}{\partial\,\epsilon_{\gamma}\,\partial\,\epsilon_{\lambda}}\bigg{|}_{\bm{\epsilon}=0}.$
(3)

This definition is only valid for the description of an initially unstressed
crystal. For a crystal under arbitrary constant stress, it is necessary to
introduce a stiffness tensor under load, $B_{ijkl}$. As hydrostatic pressure
does not reduce the crystal symmetry, it has the same symmetry as the
stiffness tensor $\bm{C}$, and can be represented in Voigt notation,
henceforth denoted as $\bm{B}$. The appropriate thermodynamical potential to
describe the stressed system is the enthalpy $H=E+p_{0}V$, where $p_{0}$
represents the initial constant pressure. The stiffness tensor under external
hydrostatic pressure can then be calculated as

$B_{\gamma\lambda}=\frac{1}{\widetilde{V}_{0}}\,\frac{\partial^{2}\,H(\bm{\epsilon})}{\partial\,\epsilon_{\gamma}\,\partial\,\epsilon_{\lambda}}\bigg{|}_{\bm{\epsilon}=0}\,,$
(4)

where $\widetilde{V}_{0}$ is the volume per unit cell for the initially
stressed structure. To derive the relation between $C_{\gamma\lambda}$ and
$B_{\gamma\lambda}$, we have to transform back to Cartesian indices and start
with arbitrary constant external stress $\sigma_{ij}$. $B_{ijkl}$ can then be
expressed as [26, 23, 27]:

$\displaystyle\begin{split}B_{ijkl}&=\widetilde{C}_{ijkl}+\\\
&\frac{1}{2}(\delta_{ik}\,\sigma_{jl}+\delta_{jk}\,\sigma_{il}+\delta_{il}\,\sigma_{jk}+\delta_{jl}\,\sigma_{ik}-2\,\delta_{kl}\,\sigma_{ij}).\end{split}$
(5)

Note that now one has to evaluate the stiffness tensor $\widetilde{C}_{ijkl}$
for the initially stressed crystal configuration at volume
$\widetilde{V}_{0}$. In this work, we consider the $\beta$ phase under
external hydrostatic pressure, $\sigma_{ij}=-P\,\delta_{ij}$ ($P>0$ for
tension). Then, the stiffness tensor under load can be written as

$\displaystyle
B_{ijkl}=\widetilde{C}_{ijkl}-P(\delta_{ik}\,\delta_{jl}+\delta_{jk}\,\delta_{il}-\delta_{kl}\,\delta_{ij}).$
(6)

Both tensors have the same symmetry, and we can thus go back to Voigt
notation. It follows that

$\displaystyle
B_{\gamma\lambda}=\tilde{C}_{\gamma\lambda}+\begin{pmatrix}-P&\phantom{-}P&\phantom{-}P&\phantom{-}0&\phantom{-}0&\phantom{i}0\\\
\phantom{-}P&-P&\phantom{-}P&\phantom{-}0&\phantom{-}0&\phantom{i}0\\\
\phantom{-}P&\phantom{-}P&-P&\phantom{-}0&\phantom{-}0&\phantom{i}0\\\
\phantom{-}0&\phantom{-}0&\phantom{-}0&-P&\phantom{-}0&\phantom{i}0\\\
\phantom{-}0&\phantom{-}0&\phantom{-}0&\phantom{-}0&-P&\phantom{i}0\\\
\phantom{-}0&\phantom{-}0&\phantom{-}0&\phantom{-}0&\phantom{-}0&-P\;\;\;\end{pmatrix}.$
(7)

The Born stability criterion [28] enables us to determine the elastic
stability of an unstressed crystal: A crystal with arbitrary symmetry is
stable if the stiffness tensor $\bm{C}$ is positive definite [28, 29]. This
condition is equivalent to $\bm{C}$ being a symmetric tensor having only
positive eigenvalues that can be calculated with standard algebraic
techniques. While the Born stability criterion is formulated for an unstressed
crystal, they can also be generalized to the case of constant external
hydrostatic load; in this case, the stiffness tensor under load, $\bm{B}$, has
to be positive definite [29].

##  III Computational details

All calculations are performed in the framework of density-functional theory
using the all-electron full-potential code exciting [30], which applies the
linearized augmented planewave plus local orbital method. Exchange-correlation
effects are treated within the generalized gradient approximation,
specifically the functional PBEsol [31] which shows high accuracy in
determining elastic properties for solids [32]. Total energies are calculated
using an 8$\times$8$\times$8 (6$\times$6$\times$6) k-grid and a planewave cut-
off $R_{\text{MT}}^{\text{min}}\,G_{\text{max}}=9.0$ for $\alpha$-Ga2O3
($\beta$-Ga2O3). We employ muffin-tin radii of
${R_{\text{MT}}^{\text{Ga}}=$1.65\text{\,}\mathrm{\text{ $a_{0}$
}}$}~{}($1.75\text{\,}\mathrm{\text{ $a_{0}$ }}$)$ and
${R_{\text{MT}}^{\text{O}}=$1.45\text{\,}\mathrm{\text{ $a_{0}$ }}$}$ for
gallium and oxygen, respectively, in $\alpha$-Ga2O3 ($\beta$-Ga2O3), where
$\mathrm{\text{ $a_{0}$ }}$ is the Bohr radius. The internal atomic positions
are relaxed until the atomic forces are smaller than 0.2 mHa $a_{0}^{-1}$.
These parameters ensure a numerical precision of
${10}^{-2}\text{\,}\mathrm{\text{ $a_{0}$ }}^{3}$ for the equilibrium volume
$V_{0}$, ${10}^{-2}\text{\,}\mathrm{GPa}$ for the bulk modulus $B_{0}$, and
${10}^{-3}$ for its pressure derivative $B^{\prime}_{0}$.

In order to calculate the full stiffness tensor at zero pressure, a total of
13 (6) different deformation types are applied to $\beta$-Ga2O3
($\alpha$-Ga2O3). The reader is referred to Ref. 25 for a full list of the
employed deformation types. For each of them, several equally spaced strain
points around the origin are created for physical strain up to
$4.5\text{\,}\mathrm{\char 37\relax}$. For every deformed structure, the
internal atomic positions are relaxed until the atomic forces are smaller than
0.2 mHa $a_{0}^{-1}$. The preparation of the deformed structures and
evaluation of second-order derivatives shown in Eqs. (3) and (4), are
performed using the ElaStic tool [25]. We create eight deformed $\beta$-Ga2O3
structures to analyze the variation of the elastic constants under strain. The
corresponding pressure values are extracted from the energy-vs-volume fit and
range from $-10$ to $35\text{\,}\mathrm{GPa}$.

All input and output files are available at the NOMAD Repository [33, 34] with
the following DOI https://dx.doi.org/10.17172/NOMAD/2021.12.02-1.

##  IV Results and Discussion

###  IV.1 Structural properties

Table 1: Calculated lattice parameters, $a$, $b$, and $c$, (in
$\mathrm{\SIUnitSymbolAngstrom}$), monoclinic angle $\beta$ (in
$\mathrm{\SIUnitSymbolDegree}$), bulk moduli, $B_{0}$ (in $\mathrm{GPa}$), and
their pressure derivatives, $B^{\prime}_{0}$, of $\alpha$-Ga2O3 and
$\beta$-Ga2O3 at zero pressure.

$\alpha$-Ga2O3 | $a$ | $c$ | ${B}_{0}$ |  $B^{\prime}_{0}$
---|---|---|---|---
Present work | 5.01 | 13.47 | 218 | 4.5
Theory [19] (LDA) | 4.95 | 13.32 | 244 | 3.8
Theory [20] (AM05) | 5.00 | 13.45 | 215 | 4.5
Experiment [14] | 4.98 | 13.43 | 252 | 4

$\beta$-Ga2O3 | $a$ | $b$ | $c$ | $\beta$ | ${B}_{0}$ |  $B^{\prime}_{0}$
---|---|---|---|---|---|---
Present work | 12.30 | 3.05 | 5.82 | 103.7 | 169 | 3.9
Theory [35] (LDA) | 12.21 | 3.03 | 5.75 | 103.6 | 219 | 3.2
Theory [20] (AM05) | 12.30 | 3.05 | 5.81 | 103.7 | 165 | 3.8
Experiment [14] | 12.23 | 3.04 | 5.80 | 103.8 | 184 | 4

The lattice parameters, bulk moduli, $B_{0}$, and their pressure derivatives,
$B^{\prime}_{0}$, are obtained as fitting parameters from the Birch-Murnaghan
equation of state (EOS). The results for the $\alpha$ and $\beta$ phase are
given in Table 1. All parameters, obtained at zero temperature and pressure,
show excellent agreement with previously published theoretical and
experimental work. The bulk modulus of the $\alpha$ phase is larger than that
of the $\beta$ phase by about $50\text{\,}\mathrm{GPa}$. Considering that
$\alpha$-Ga2O3 is the more compressed phase and the bulk modulus is a measure
of resistance against uniform compression, such a result is expected. Compared
to experiment, it is underestimated by about 10 - 20 $\mathrm{\char
37\relax}$, which is within the typical accuracy of a semi-local DFT
calculation of elastic properties [36]. Note that $B_{0}$ here is obtained
from a fit and not explicitly calculated as a linear combination of second-
order elastic constants (see Section IV.3).

###  IV.2 Phase transition

Using the energy-vs-volume relation at high pressures, we can further analyze
the phase transition from the $\beta$ to the $\alpha$ phase. The Gibbs free
energy, $G=U+pV-TS$, dictates the structural stability at a given temperature
and pressure. The complete calculation of this quantity would require the full
phonon spectrum. In this work, we focus on the enthalpy $H=E+pV$, i.e., the
low-temperature case, where the internal energy, $U\approx E$, is determined
by the Birch-Murnaghan EOS. Such an approach is justified for the pressure-
induced properties of hard materials [37]. For a given pressure, the crystal
phase with the lowest enthalpy is the most stable one while a crossing point
between two phases indicates a first-order phase transition. This transition
pressure is purely obtained by thermodynamical considerations and will
henceforth be denoted as $p_{\text{t}}$.

Figure 2: Equations of state for the $\alpha$ and $\beta$ phase of Ga2O3.
Energies are relative to values of $\beta$-Ga2O3. The dots and crosses are the
computed data points for $\beta$-Ga2O3 and $\alpha$-Ga2O3, respectively. The
solid and dashed lines represent the corresponding Birch-Murnaghan (B-M) fits.
The dotted line represents the common tangent of the EOS for both phases. The
transition pressure $p_{\text{t}}$ is determined by the common gradient.

The calculated Birch-Murnaghan EOS for both phases are given in Fig. 2. We
obtain a transition pressure of $p_{\text{t}}=$2.6\text{\,}\mathrm{GPa}$$.
This value can be compared with the experimental transition onset, owing to
the fact that a new phase is thermodynamically favored once this pressure is
reached. We therefore consider $p_{\text{t}}$ as the lower bound for a phase
transition setting in. As illustrated in the Introduction, the discrepancy
between published transition pressures is quite high as the experimental
conditions play a vital role in determining the transition onset, where
temperature, pressure medium, sample size and type can have a dramatic impact.
Consequently, reported experimental values range from 3 to
$20\text{\,}\mathrm{GPa}$ [14, 16, 38, 17, 13]. However, also theoretical
values range from $2$ to $17\text{\,}\mathrm{GPa}$ [18, 19, 20, 21]. Here,
differences can be assigned to the usage of different exchange-correlation
functionals as well as as the treatment of the Ga 3$d$ states (as core or
valence states) in pseudopotential approaches.

###  IV.3 Elastic stability at different pressures

####  IV.3.1 Ambient pressure

Table 2: Calculated second-order elastic constants at ambient pressure for the
$\beta$ phase of Ga2O3, given in units of $\mathrm{GPa}$, compared to
previously published results.

$\beta$-Ga2O3 | $C_{11}$ | $C_{12}$ | $C_{13}$ | $C_{15}$ | $C_{22}$ | $C_{23}$ | $C_{25}$ | $C_{33}$ | $C_{35}$ | $C_{44}$ | $C_{46}$ | $C_{55}$ |  $C_{66}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Present work | 220.5 | 116.4 | 128.3 | -17.3 | 329.5 | 79.1 | 11.4 | 326.8 | 7.2 | 50.0 | 18.1 | 66.1 |  191.5
Theory [20] | 223.1 | 116.5 | 125.3 | -17.4 | 333.2 | 75.0 | 12.2 | 330.0 | 7.3 | 50.3 | 17.4 | 68.6 |  194.2
Expt. [39] | 242.8 | 128.0 | 160.0 | -11.6 | 343.8 | 70.9 |  10.4 | 347.4 | 1.0 | 47.8 |  15.6 | 88.6 | 104.0

Table 3: Calculated second-order elastic constants at ambient pressure for the
$\alpha$ phase of Ga2O3, given in units of $\mathrm{GPa}$, compared to
previously published theoretical results.

$\alpha$-Ga2O3 | $C_{11}$ | $C_{12}$ | $C_{13}$ | $C_{14}$ | $C_{33}$ |  $C_{44}$
---|---|---|---|---|---|---
Present work | 380.3 | 174.5 | 128.5 | -16.6 | 342.9 | 80.1
Theory [20] | 381.5 | 173.6 | 126.0 | -17.3 | 345.8 | 79.7

Depending on the crystal symmetry, the stiffness tensor can have up to 21
independent components. In the case of the monoclinic $\beta$ phase
(rhombohedral $\alpha$ phase), this number reduces to 13 (6). The calculated
second-order elastic constants at zero pressure for both phases are summarized
in Tables 3 and 3, showing excellent agreement with other theoretical results.
To compare our results to those of Ref. 20, we employ the same lattice
representation. For $\beta$-Ga2O3 it is chosen such that the $y$ axis is
parallel to $b$, the $x$ axis is parallel to $a$, and the $c$ axis lies in the
$x$-$z$ plane (see Fig. 1).

Table 4: Calculated elastic moduli (in units of $\mathrm{GPa}$) given as Voigt’s, Reuss’s, and Hill’s average. For the bulk modulus we also show the value obtained from the Birch-Murnaghan EOS in Fig. (2). $\beta$-Ga2O3 | ${B}$ | ${G}$ | ${E}$ |  ${\nu}$
---|---|---|---|---
Voigt | 169.4 | 78.4 | 203.7 | 0.30
Reuss | 166.7 | 66.2 | 175.3 | 0.32
Hill | 168.0 | 72.3 | 189.6 | 0.31
Birch-Murnaghan EOS | 168.8 |  |  |
$\alpha$-Ga2O3 | ${B}$ | ${G}$ | ${E}$ |  ${\nu}$
Voigt | 218.5 | 97.4 | 254.5 | 0.31
Reuss | 216.4 | 92.2 | 242.0 | 0.31
Hill | 217.4 | 94.8 | 248.2 | 0.31
B-M EOS | 218.4 |  |  |

We also calculate the elastic moduli as linear combinations of second-order
elastic constants and compare the bulk modulus with the one obtained from the
Birch-Murnaghan EOS, $B_{0}$, shown in Table 4. The values for the bulk moduli
from both methods show excellent agreement, further validating the precision
of our calculations. All elastic moduli of the $\beta$ phase are smaller than
those of the $\alpha$ phase. This is expected, since the latter is obtained
from compressing $\beta$-Ga2O3 under hydrostatic strain, and a higher density
leads to more resistance against strain and, therefore, larger elastic moduli.

The second-order elastic constants of both phases exhibit pronounced
anisotropy due to crystal symmetry. The diagonal terms, $C_{11}$, $C_{22}$,
and $C_{33}$, have the highest values in excess of $200\text{\,}\mathrm{GPa}$.
As large diagonal components mean a high degree of hardness against strain in
the respective directions, both phases strongly resist deformations along the
main axes. This situation is reversed for the shear-strain components (indices
$4$ to $6$). The calculated values are much smaller than the axial strain
components (indices $1$ to $3$). They are as low as
$C_{25}=$11.4\text{\,}\mathrm{GPa}$$ and $C_{35}=$7.2\text{\,}\mathrm{GPa}$$
for the $\beta$ phase, and $C_{14}=$-16.5\text{\,}\mathrm{GPa}$$ for the
$\alpha$ phase. These findings suggest that both phases are susceptible to
shear strains. The elastic moduli further validate this assumption, as the
shear modulus for both phases is smaller than the bulk and Young modulus,
i.e., $G<B<E$.

Importantly, for both the $\beta$ and $\alpha$ phase, all eigenvalues of the
stiffness tensors are positive. Thus, according to the Born stability
criterion, they are elastically stable at equilibrium at
$0\text{\,}\mathrm{K}$. This coincides with previously published results [11],
where $\beta$-Ga2O3 was identified as the thermodynamically stable phase and
$\alpha$-Ga2O3 as a metastable phase.

####  IV.3.2 Variation under hydrostatic pressure

We now explore how the second-order elastic constants of the $\beta$ phase
react to hydrostatic pressure. To this extent, we perform calculations on
configurations that are obtained by isotropically straining the crystal, by
applying the transformation

$\bm{R}(\epsilon)=(1+\epsilon)\bm{R}_{P=0}\,,$

where $\epsilon$ is a constant value, and calculate their stiffness tensors
under load, $\bm{B}$ (Eq. (4)). The chosen values for $\epsilon$ correspond to
pressure values between $-10\text{\,}\mathrm{GPa}$ and
$35\text{\,}\mathrm{GPa}$. The results are summarized in Table 5. Overall, the
elastic constants increase in value with increasing pressure, reflecting that
the denser structures are more resistant to strain. The eigenvalues of the
stiffness tensor are positive up to $20\text{\,}\mathrm{GPa}$, indicating
mechanical stability of the $\beta$ phase under strain. The structure at
$35\text{\,}\mathrm{GPa}$ is the only unstable one. We conclude that the
critical pressure $p_{\text{c}}$ must thus be much higher than the
thermodynamical transition pressure of $2.6\text{\,}\mathrm{GPa}$ obtained in
Section IV.2. As such, $p_{\text{c}}$ can be seen as the upper bound of the
phase transition pressure.

Table 5: Second-order elastic constants of $\beta$-Ga2O3 as a function of hydrostatic pressure $P$. Pressure and elastic constants are given in units of $\mathrm{GPa}$. $V_{0}$ is the volume of the unstressed crystal. $P$ | $V/V_{0}$ | $B_{11}$ | $B_{12}$ | $B_{13}$ | $B_{15}$ | $B_{22}$ | $B_{23}$ | $B_{25}$ | $B_{33}$ | $B_{35}$ | $B_{44}$ | $B_{46}$ | $B_{55}$ |  $B_{66}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
-10.0 | 1.07 | 175.9 |  190.6 |  184.9 | -15.0 | 272.3 | 37.8 | 13.8 | 277.3 | 5.3 | 25.8 |  17.6 | 66.7 | 87.4
-5.0 | 1.03 | 201.5 | 103.7 | 106.0 | -15.2 | 302.4 | 58.2 | 12.9 | 303.0 | 6.6 | 40.9 | 13.2 | 66.7 | 90.0
0.0 | 1.00 | 220.5 | 116.4 | 128.3 | -17.3 | 329.5 | 79.1 | 11.4 | 326.8 | 7.2 | 50.0 | 18.1 | 66.1 | 91.5
0.5 | 0.99 | 222.2 | 117.2 | 131.7 | -18.4 | 332.2 | 79.6 | 10.3 | 328.8 | 7.2 | 50.9 | 18.2 | 65.8 | 91.7
2.5 | 0.99 | 229.3 | 122.6 | 140.0 | -18.8 | 343.4 | 87.7 | 10.8 | 336.3 | 7.7 | 53.6 | 19.4 | 65.4 | 92.1
4.9 | 0.97 | 236.3 | 129.9 | 149.7 | -19.0 | 357.1 | 98.5 | 11.5 | 343.9 | 8.8 | 57.1 | 20.2 | 64.9 | 92.5
10.0 | 0.95 | 249.8 | 143.7 | 174.5 | -23.9 | 382.5 | 118.4 | 8.6 | 361.7 | 9.2 | 61.0 | 22.5 | 62.4 | 92.7
19.6 | 0.91 | 265.2 | 170.2 | 222.3 | -33.7 | 429.7 | 160.3 | 4.1 | 392.2 | 10.4 | 68.3 | 24.2 | 56.5 | 91.4
35.0 | 0.86 | 279.5 | 212.5 | 307.0 | -52.9 | 502.4 | 229.5 | -3.0 | 433.7 | 11.5 | 78.9 | 22.2 | 43.7 | 84.6
Figure 3:  Instability criteria $M(P)$ and $N(P)$, as defined in Eq. (9), as a
function of hydrostatic pressure $P$. The red dots represent the calculated
data points for the strained $\beta$-Ga2O3 structures. The solid lines
indicate polynomial fits of order 2 and 4 to the data points for $M(P)$ and
$N(P)$, respectively.

In order to estimate $p_{\text{c}}$, we further analyze our results using the
Born stability criterion. As an alternative to calculating the eigenvalues of
the stiffness tensor, we aim at finding a closed mathematical expression for
its conditions. [29, 40, 41]. This can be done, for example, by making use of
the leading principal minors of the stiffness tensor [29]. The onset pressure
of mechanical instabilities can be estimated by evaluating these expressions
over a certain pressure range. For cubic systems under hydrostatic stress, the
second-order constants $B_{\gamma\delta}$ only differ from the stiffness
constants $C_{\gamma\delta}$ by a linear term in $P$ [29, 40, 41]. Such linear
equations are not attainable for a monoclinic system. As an alternative
approach, we evaluate the equation

$\det\left|\bm{B}\right|=0\,,$ (8)

to explore when the system becomes unstable [41]. This enables us to estimate
the critical pressures only from the stiffness tensor under load. In the case
of $\beta$-Ga2O3, Eq. (8) is a polynomial equation of 6th degree in the
pressure which can be expressed in the following form:

$\displaystyle\underbrace{\left(P^{2}-P\,(C_{44}+C_{66})+C_{44}\,C_{66}-C_{46}^{2}\right)}_{M(P)}\cdot\underbrace{\mathcal{P}_{4}(P)}_{N(P)}$
$\displaystyle=0\,.$ (9)

The first term, denoted as $M(P)$, includes only shear components of the
stiffness tensor and can be identified as a shear instability criterion. The
second term, $\mathcal{P}_{4}(P)$, is a polynomial of order four in $P$
(denoted as $N(P)$) which is provided in the Appendix. The structure is
mechanically unstable if one term equals zero. We now estimate the critical
pressure(s) by interpolating our results for $M(P)$ and $N(P)$. By doing so,
we can identify from the zero of $M(P)$ whether the transition is solely
occurring due to shear strain. Our interpolation of both expressions is to be
used with caution for very high ($>$50\text{\,}\mathrm{GPa}$$) and very low
pressures ($<$-20\text{\,}\mathrm{GPa}$$). We want to emphasize that the aim
is to estimate the critical pressure in the experimentally relevant range of
up to $50\text{\,}\mathrm{GPa}$. The results are illustrated in Fig. 3. We
obtain two critical pressures of $p_{\text{c}}=$-21.4\text{\,}\mathrm{GPa}$$
and $p_{\text{c}}=$32.4\text{\,}\mathrm{GPa}$$. This indicates that
$\beta$-Ga2O3 is mechanically stable in the range of
$$-21.4\text{\,}\mathrm{GPa}$<P<$32.4\text{\,}\mathrm{GPa}$$. The upper bound
agrees well with the observations from several experiments as the phase
transition is fully completed above $30$ to $40\text{\,}\mathrm{GPa}$, i.e.,
only $\alpha$-Ga2O3 is remaining in the sample [13, 14, 15]. Only in Ref. 21 a
theoretical critical pressure below $30\text{\,}\mathrm{GPa}$ was reported.
Their transition pressure of $19.4\text{\,}\mathrm{GPa}$ is, however,
drastically higher than our calculated value of $2.6\text{\,}\mathrm{GPa}$.
Note that our results are obtained at $0\text{\,}\mathrm{K}$. First-principles
calculations show that the elastic constants for $\beta$-Ga2O3 are decreasing
with temperature; however, the effect is small at room temperature [42].
Therefore, we do not expect a significant change of the critical pressure at
ambient temperature. The lower bound would indicate the emergence of another
metastable phase below $-21.4\text{\,}\mathrm{GPa}$ which is due to an
instability of the pure shear criterion, $M(P)$=0. To our knowledge, no
studies have been performed with negative pressure for any of the
sesquioxides. In principle, it is possible to reach negative pressure values
on the order of a couple of $\mathrm{GPa}$ in solids [43].

While most studies report the phase transition with quasihydrostatic pressure
mediums, the phase transition also occurs under non-hydrostatic conditions
[14]. This, in turn, indicates that an additional mechanical instability may
arise from non-uniform stress, which would reduce the symmetry of the
stiffness tensor under load, $\bm{B}$. We have so far not considered this in
our analysis.

##  V Summary

We have investigated the structural and elastic properties of Ga2O3 in the
rhombohedral $\alpha$ and monoclinic $\beta$ phase from first principles.
Based on our results, a phase transition from $\beta$\- to $\alpha$-Ga2O3 is
energetically favored at $p_{\text{t}}=$2.6\text{\,}\mathrm{GPa}$$. The
calculated full stiffness tensors, $\bm{C}$, of both phases at ambient
pressure show pronounced anisotropy as well as susceptibility to shear strain,
indicated by the small shear moduli of
$G_{V}^{\beta}=$78\text{\,}\mathrm{GPa}$$ and
$G_{V}^{\alpha}=$97\text{\,}\mathrm{GPa}$$, respectively. Investigating the
variation of the stiffness tensor under hydrostatic pressure, we observe that,
according to the Born stability criterion, $\beta$-Ga2O3 is becoming
mechanically unstable at a critical pressure of
$p_{\text{c}}=$32.4\text{\,}\mathrm{GPa}$$. The transition pressure
$p_{\text{t}}$ can be seen as a lower bound for the phase transition and
agrees well with the transition onset in previously reported experimental and
theoretical results [14, 20, 16, 18]. Only considering this pressure value
obtained from thermodynamics is, however, not sufficient to explain the full
range of experimentally observed transition pressures. While a phase
transition is energetically favored above $p_{\text{t}}$, there are additional
kinetic barriers that must be overcome. Only when the $\beta$ phase is
mechanically unstable, i.e. at $p_{\text{c}}=$32.4\text{\,}\mathrm{GPa}$$, we
expect the phase transition to be completed. This agrees well with
experimental observations showing completion of the transition only above
$30\text{\,}\mathrm{GPa}$ [14, 15, 16]. In addition, we find a critical
pressure of $-21.4\text{\,}\mathrm{GPa}$. Phonon calculations in this pressure
range could provide insight into the emergence of a novel metastable phase for
negative pressures. Experiments with negative pressures, at values below a few
$\mathrm{GPa}$ have so far not been conducted but could point to novel
discoveries in the future.

Our first-principles approach, being successfully demonstrated here for the
wide-gap oxide Ga2O3, though requiring hydrostatic pressure, is independent of
crystal symmetry, and thus could be applied to other materials. For example, a
similar analysis could thus be conducted for other sesquioxides to estimate
stability windows and the possible emergence of novel metastable phases.

##  VI Acknowledgements

This work was performed in the framework of GraFOx, a Leibniz-Science Campus
partially funded by the Leibniz Association. The authors acknowledge the
North-German Supercomputing Alliance (HLRN) for providing HPC resources that
have contributed to the research results reported in this paper. The authors
would like to thank Maria Troppenz for critically reading the manuscript.

##  VII Appendix

The term $N(P)$ in Eq. 9 is a polynomial of order 4 in $P$ and can be
expressed in terms of the elastic constants $C_{\gamma\lambda}$ as follows:

$N(P)=\hskip 4.0ptb_{4}\,P^{4}+b_{3}\,P^{3}+b_{2}\,P^{2}+b_{1}\,P+b_{0}\,,\\\
$

with

$\displaystyle b_{4}=$ $\displaystyle-4\,,$ $\displaystyle b_{3}=$
$\displaystyle-4C_{12}-4C_{13}-4C_{23}+4C_{55}\,,$ $\displaystyle b_{2}=$
$\displaystyle\hskip
4.0ptC_{11}C_{22}+2C_{11}C_{23}+C_{11}C_{33}-C_{12}^{2}-2C_{12}C_{13}-2C_{12}C_{23}+2C_{12}C_{33}+4C_{12}C_{55}-C_{13}^{2}+2C_{13}C_{22}$
$\displaystyle-2C_{13}C_{23}+4C_{13}C_{55}-4C_{15}C_{25}-4C_{15}C_{35}+C_{22}C_{33}-C_{23}^{2}+4C_{23}C_{55}-4C_{25}C_{35}\,,$
$\displaystyle b_{1}=$ $\displaystyle-
C_{11}C_{22}C_{33}-C_{11}C_{22}C_{55}+C_{11}C_{23}^{2}-2C_{11}C_{23}C_{55}+C_{11}C_{25}^{2}+2C_{11}C_{25}C_{35}-C_{11}C_{33}C_{55}$
$\displaystyle+C_{11}C_{35}^{2}+C_{12}^{2}C_{33}+C_{12}^{2}C_{55}-2C_{12}C_{13}C_{23}+2C_{12}C_{13}C_{55}-2C_{12}C_{15}C_{25}-2C_{12}C_{15}C_{35}$
$\displaystyle+2C_{12}C_{23}C_{55}-2C_{12}C_{25}C_{35}-2C_{12}C_{33}C_{55}+2C_{12}C_{35}^{2}+C_{13}^{2}C_{22}+C_{13}^{2}C_{55}-2C_{13}C_{15}C_{25}$
$\displaystyle-2C_{13}C_{15}C_{35}-2C_{13}C_{22}C_{55}+2C_{13}C_{23}C_{55}+2C_{13}C_{25}^{2}-2C_{13}C_{25}C_{35}+C_{15}^{2}C_{22}+2C_{15}^{2}C_{23}$
$\displaystyle+C_{15}^{2}C_{33}+2C_{15}C_{22}C_{35}-2C_{15}C_{23}C_{25}-2C_{15}C_{23}C_{35}+2C_{15}C_{25}C_{33}-C_{22}C_{33}C_{55}$
$\displaystyle+C_{22}C_{35}^{2}+C_{23}^{2}C_{55}-2C_{23}C_{25}C_{35}+C_{25}^{2}C_{33}\,,$
$\displaystyle b_{0}=$ $\displaystyle\hskip
4.0ptC_{11}C_{22}C_{33}C_{55}-C_{11}C_{22}C_{35}^{2}-C_{11}C_{23}^{2}C_{55}+2C_{11}C_{23}C_{25}C_{35}-C_{11}C_{25}^{2}C_{33}-C_{12}^{2}C_{33}C_{55}$
$\displaystyle+C_{12}^{2}C_{35}^{2}+2C_{12}C_{13}C_{23}C_{55}-2C_{12}C_{13}C_{25}C_{35}-2C_{12}C_{15}C_{23}C_{35}+2C_{12}C_{15}C_{25}C_{33}$
$\displaystyle-
C_{13}^{2}C_{22}C_{55}+C_{13}^{2}C_{25}^{2}+2C_{13}C_{15}C_{22}C_{35}-2C_{13}C_{15}C_{23}C_{25}-C_{15}^{2}C_{22}C_{33}+C_{15}^{2}C_{23}^{2}\,.$

## References

  * Huang  _et al._ [2006] Y. Huang, S. Yue, Z. Wang, Q. Wang, C. Shi, Z. Xu, X. D. Bai, C. Tang, and C. Gu, Preparation and electrical properties of ultrafine Ga2O3 nanowires, The Journal of Physical Chemistry B 110, 796 (2006).
  * Arnold  _et al._ [2009] S. P. Arnold, S. M. Prokes, F. K. Perkins, and M. E. Zaghloul, Design and performance of a simple, room-temperature Ga2O3 nanowire gas sensor, Applied Physics Letters 95, 103102 (2009).
  * Cuong  _et al._ [2009] N. D. Cuong, Y. W. Park, and S. G. Yoon, Microstructural and electrical properties of Ga2O3 nanowires grown at various temperatures by vapour-liquid-solid technique, Sensors and Actuators B: Chemical 140, 240 (2009).
  * Liu  _et al._ [2008] Z. Liu, T. Yamazaki, Y. Shen, T. Kikuta, N. Nakatani, and Y. Li, O2 and CO sensing of Ga2O3 multiple nanowire gas sensors, Sensors and Actuators B: Chemical 129, 666 (2008).
  * Mazeina  _et al._ [2010] L. Mazeina, F. K. Perkins, V. M. Bermudez, S. P. Arnold, and S. M. Prokes, Functionalized Ga2O3 nanowires as active material in room temperature capacitance-based gas sensors, Langmuir 26, 13722 (2010).
  * Chang  _et al._ [2005] P.-C. Chang, Z. Fan, W.-Y. Tseng, A. Rajagopal, and J. G. Lu, $\beta$-Ga2O3 nanowires: Synthesis, characterization, and p-channel field-effect transistor, Applied Physics Letters 87, 222102 (2005).
  * Feng  _et al._ [2006] P. Feng, J. Y. Zhang, Q. H. Li, and T. H. Wang, Individual $\beta$-Ga2O3 nanowires as solar-blind photodetectors, Applied Physics Letters 88, 153107 (2006).
  * Li  _et al._ [2011] L. Li, E. Auer, M. Liao, X. Fang, T. Zhai, U. K. Gautam, A. Lugstein, Y. Koide, Y. Bando, and D. Golberg, Deep-ultraviolet solar-blind photoconductivity of individual gallium oxide nanobelts, Nanoscale 3, 1120 (2011).
  * Kokubun  _et al._ [2007] Y. Kokubun, K. Miura, F. Endo, and S. Nakagomi, Sol-gel prepared $\beta$-Ga2O3 thin films for ultraviolet photodetectors, Applied Physics Letters 90, 031912 (2007).
  * Weng  _et al._ [2010] W. Y. Weng, T. J. Hsueh, S. J. Chang, G. J. Huang, and S. P. Chang, A solar-blind $\beta$-Ga2O3 nanowire photodetector, IEEE Photonics Technology Letters 22, 709 (2010).
  * Roy  _et al._ [1952] R. Roy, V. G. Hill, and E. F. Osborn, Polymorphism of Ga2O3 and the system Ga2O3$-$H2O, Journal of the American Chemical Society 74, 719 (1952).
  * Remeika and Marezio [1966] J. P. Remeika and M. Marezio, Growth of $\alpha$‐Ga2O3 single crystals at 44 kbars, Applied Physics Letters 8, 87 (1966).
  * Lipinska-Kalita  _et al._ [2003] K. E. Lipinska-Kalita, B. Chen, M. B. Kruger, Y. Ohki, J. Murowchick, and E. P. Gogol, High-pressure x-ray diffraction studies of the nanostructured transparent vitroceramic medium K2O$-$SiO2$-$Ga2O3, Phys. Rev. B 68, 035209 (2003).
  * Lipinska-Kalita  _et al._ [2008] K. E. Lipinska-Kalita, P. E. Kalita, O. A. Hemmers, and T. Hartmann, Equation of state of gallium oxide to $70\phantom{\rule{3.00003pt}{0.0pt}}\mathrm{GPa}$: Comparison of quasihydrostatic and nonhydrostatic compression, Phys. Rev. B 77, 094123 (2008).
  * Wang  _et al._ [2010] H. Wang, Y. He, W. Chen, Y. W. Zeng, K. Stahl, T. Kikegawa, and J. Z. Jiang, High-pressure behavior of $\beta$-Ga2O3 nanocrystals, Journal of Applied Physics 107, 033520 (2010).
  * Machon  _et al._ [2006] D. Machon, P. F. McMillan, B. Xu, and J. Dong, High-pressure study of the $\beta$-to-$\alpha$ transition in Ga2O3, Phys. Rev. B 73, 094125 (2006).
  * Kishimura and Matsumoto [2018] H. Kishimura and H. Matsumoto, Evaluation of the shock-induced phase transition in $\beta$-Ga2O3, Japanese Journal of Applied Physics 57, 125503 (2018).
  * Kroll [2005] P. Kroll, Spinel-type gallium oxynitrides attainable at high pressure and high temperature, Phys. Rev. B 72, 144407 (2005).
  * Litimein  _et al._ [2009] F. Litimein, D. Rached, R. Khenata, and H. Baltache, FPLAPW study of the structural, electronic, and optical properties of Ga2O3: Monoclinic and hexagonal phases, Journal of Alloys and Compounds 488, 148 (2009).
  * Furthmüller and Bechstedt [2016] J. Furthmüller and F. Bechstedt, Quasiparticle bands and spectra of Ga2O3 polymorphs, Phys. Rev. B 93, 115204 (2016).
  * Luan  _et al._ [2019] S. Luan, L. Dong, and R. Jia, Analysis of the structural, anisotropic elastic and electronic properties of $\beta$-Ga2O3 with various pressures, Journal of Crystal Growth 505, 74 (2019).
  * Timoshenko [1951] S. Timoshenko,  _Theory of Elasticity_ (McGraw-Hill, 1951).
  * Wallace [1972] D. C. Wallace,  _Thermodynamics of Crystals_ (Wiley, 1972).
  * Landau and Lifschitz [1986] L. D. Landau and E. M. Lifschitz,  _Theory of Elasticity_ (Elsevier, 1986).
  * Golesorkhtabar  _et al._ [2013] R. Golesorkhtabar, P. Pavone, J. Spitaler, P. Puschnig, and C. Draxl, Elastic: A tool for calculating second-order elastic constants from first principles, Comp. Phys. Commun. 184, 1861 (2013).
  * Grimvall  _et al._ [2012] G. Grimvall, B. Magyari-Köpe, V. Ozoliņš, and K. A. Persson, Lattice instabilities in metallic elements, Rev. Mod. Phys. 84, 945 (2012).
  * Wallace [1967] D. C. Wallace, Thermoelasticity of stressed materials and comparison of various elastic constants, Phys. Rev. 162, 776 (1967).
  * Born [1940] M. Born, On the stability of crystal lattices. I, Math. Proc. Cambridge Philos. Soc. 36, 160 (1940).
  * Mouhat and Coudert [2014] F. Mouhat and F.-X. Coudert, Necessary and sufficient elastic stability conditions in various crystal systems, Phys. Rev. B 90, 224104 (2014).
  * Gulans  _et al._ [2014] A. Gulans, S. Kontur, C. Meisenbichler, D. Nabok, P. Pavone, S. Rigamonti, S. Sagmeister, U. Werner, and C. Draxl, exciting: a full-potential all-electron package implementing density-functional theory and many-body perturbation theory, J.Phys.: Condens. Matter 26, 363202 (2014).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the density-gradient expansion for exchange in solids and surfaces, Phys. Rev. Lett. 100, 136406 (2008).
  * Râsander and Moram [2015] M. Râsander and M. A. Moram, On the accuracy of commonly used density functional approximations in determining the elastic constants of insulators and semiconductors, J. Chem. Phys. 143, 144104 (2015).
  * Draxl and Scheffler [2018] C. Draxl and M. Scheffler, NOMAD: The FAIR concept for big-data-driven materials science, MRS Bulletin 43, 676 (2018).
  * Draxl and Scheffler [2019] C. Draxl and M. Scheffler, The NOMAD laboratory: From data sharing to artificial intelligence, J. Phys. Mater. 2, 036001 (2019).
  * He  _et al._ [2006] H. He, R. Orlando, M. A. Blanco, R. Pandey, E. Amzallag, I. Baraille, and M. Rérat, First-principles study of the structural, electronic, and optical properties of Ga2O3 in its monoclinic and hexagonal phases, Phys. Rev. B 74, 195123 (2006).
  * Giustino [2014] F. Giustino,  _Materials Modelling using Density Functional Theory: Properties and Predictions_ (Oxford University Press, 2014).
  * Karch  _et al._ [1996] K. Karch, F. Bechstedt, P. Pavone, and D. Strauch, Pressure-dependent properties of SiC polytypes, Phys. Rev. B 53, 13400 (1996).
  * Yan-Mei  _et al._ [2008] M. Yan-Mei, C. Hai-Yong, Y. Kai-Feng, L. Min, C. Qi-Liang, L. Jing, and Z. Guang-Tian, High-pressure and high-temperature behaviour of gallium oxide, Chin. Phys. Lett. 25, 1603 (2008).
  * Adachi  _et al._ [2018] K. Adachi, H. Ogi, N. Takeuchi, N. Nakamura, H. Watanabe, T. Ito, and Y. Ozaki, Unusual elasticity of monoclinic $\beta$-Ga2O3, Journal of Applied Physics 124, 085102 (2018).
  * Wang  _et al._ [1993] J. Wang, S. Yip, S. R. Phillpot, and D. Wolf, Crystal instabilities at finite strain, Phys. Rev. Lett. 71, 4182 (1993).
  * Wang  _et al._ [1995] J. Wang, J. Li, S. Yip, S. Phillpot, and D. Wolf, Mechanical instabilities of homogeneous crystals, Phys. Rev. B 52, 12627 (1995).
  * Santia  _et al._ [2019] M. Santia, N. Tandon, and J. Albrecht, Effects of strain on the optical and thermodynamic properties of $\beta$-Ga2O3, Solid State Communications 297, 1 (2019).
  * Imre [2007] A. R. Imre, On the existence of negative pressure states, physica status solidi (b) 244, 893 (2007).