# Microwave graphs analogs for the voltage drop in three-terminal devices with
orthogonal, unitary and symplectic symmetry

F. Castañeda-Ramírez  Departamento de Física, Universidad Autónoma
Metropolitana-Iztapalapa, Apartado Postal 55-534, 09340 Ciudad de México,
Mexico  A. M. Martínez-Argüello<EMAIL_ADDRESS>Instituto de
Ciencias Físicas, Universidad Nacional Autónoma de México, Apartado Postal
48-3, 62210 Cuernavaca, Mor., Mexico  T. Hofmann  Fachbereich Physik der
Philipps-Universität Marburg, D-35032 Marburg, Germany  A. Rehemanjiang
Fachbereich Physik der Philipps-Universität Marburg, D-35032 Marburg, Germany
M. Martínez-Mares  Departamento de Física, Universidad Autónoma Metropolitana-
Iztapalapa, Apartado Postal 55-534, 09340 Ciudad de México, Mexico  J. A.
Méndez-Bermúdez  Instituto de Física, Benemérita Universidad Autónoma de
Puebla, Apartado Postal J-48, 72570 Puebla, Pue., Mexico  U. Kuhl  Université
Côte d’Azur, CNRS, Institut de Physique de Nice (INPHYNI), 06108 Nice, France,
EU  H.-J. Stöckmann  Fachbereich Physik der Philipps-Universität Marburg,
D-35032 Marburg, Germany

###### Abstract

Transmission measurements through three-port microwave graphs are performed in
a symmetric setting, in analogy to three-terminal voltage drop devices with
orthogonal, unitary, and symplectic symmetry. The terminal used as a probe is
symmetrically located between two chaotic graphs, each graph is connected to
one port, the input and the output, respectively. The analysis of the
experimental data exhibit the weak localization and antilocalization phenomena
in a clear fashion. We find a good agreement with theoretical predictions,
provided that the effect of dissipation and imperfect coupling to the ports
are taken into account.

###### pacs:

73.23.-b, 73.21.Hb, 72.10.-d, 72.15.Rn

##  I Introduction

The wave origin of quantum interference in many physical phenomena opens the
possibility that analogous classical wave systems may emulate quantum devices
Fyodorov2005a . For example, classical wave systems have been used as
auxiliary tools to understand transport properties of multiterminal quantum
systems. The dissipation that occurs in classical wave systems and the
imperfect coupling to the leads that feed the system have not been
disadvantageous, but interesting phenomena deserving to be studied to analyze
their effects in transport properties Schanze2001 ; Schanze2005 ; Schafer2003
; Kuhl2005 .

An opportunity to study many electrical or thermal conduction properties of
solid state physics by classical wave systems has been opened once the problem
of electrical or thermal conduction is reduced to a scattering problem. For
instance, Landauer formula states that the electrical conductance is
proportional to the transmission coefficient Buttiker1988 ; an equivalent
Landauer formula is valid for thermal conductance Schwab2000 . Therefore, a
lot of research has been devoted to the study of conduction through two-
terminal configurations from the theoretical and experimental points of view,
using quantum mechanics Buttiker1984 ; Brouwer1997 ; Keller1996 ; Marcus1992 ;
Chan1995  as well as classical wave physics Schanze2001 ; Schanze2005 ;
Moises2005 ; Enrique2016 .

Multiterminal systems have also been studied in Buttiker1986 ; Gopar1994 ;
Arrachea2008 ; Texier2016 ; Foieri2009 ; DAmato1990 ; Cattena2014 ; Song1998 ;
Gao2005 . Among them three-terminal systems have been considered in which the
voltage drop along the system is the observable of interest Angel2018 ;
Angel2019 . While two of the terminals are connected to fixed electronic
reservoirs to feed the system, the third one is used as a probe that tunes the
voltage drop. The voltage drop along a disordered wire was analyzed for the
one-channel case Godoy1992a ; Godoy1992b . It has been shown theoretically
that the value of voltage drop lies between those on the two terminals; in
fact, it consists of the average of the voltages in the terminals and a
deviation term, which contains specific information of the system through the
transmsission coefficients from the terminals to the probe. This _voltage drop
deviation_ term, $f$, is an important quantity in the study of the voltage
drop reducing the problem to a scattering problem. It has been shown that $f$
fluctuates from sample-to-sample and its statistical distribution has a
completely different shape for metallic and insulating regimes Godoy1992a ;
Godoy1992b . The same quantity was also considered for a chaotic system,
numerically simulated using Random Matrix Theory (RMT) Angel2014 .

More recently, an experiment with microwave graphs has been performed for an
asymmetric configuration, where the probe is located on one side of a chaotic
graph Angel2018 . Analytical and numerical procedures were carried out to
compute $f$. Remarkably, the results show significant differences with respect
to the corresponding disordered case. Here, we deepen in the understanding of
the voltage drop by proposing an experiment with microwave graphs but locating
the probe between two chaotic graphs, as shown in Fig. 1. Our proposal is
accompanied by analytical predictions.

The paper is organized as follows. In the next section we obtain the
analytical expressions for the voltage drop for a measurement at the middle of
two scattering devices. It is written in terms of the scattering matrices of
the individual devices. In section III the analytical procedure, using RMT
calculations for the scattering matrices, is developed for the statistical
distribution of the deviation $f$ assuming chaotic scattering devices. The
experimental realization with microwave graphs is explained in section IV,
where we compare experimental results with theoretical predictions. Finally,
we present our conclusions in section V.

##  II Voltage drop reduced to a scattering problem

###  II.1 Voltage drop at the middle of a symmetric device

To measure the voltage drop along a quantum wire in the simplest
configuration, a three-terminal system is needed in which one of the ports is
used as a probe by tuning its voltage to zero current. A schematic view of the
system where the probe is located between two quantum devices is shown in Fig.
1. We are interested in the one-channel situation of the leads connecting the
system, for which the zero current in the probe implies that Buttiker1988

$\mu_{3}=\frac{1}{2}\left(\mu_{1}+\mu_{2}\right)+\frac{1}{2}\left(\mu_{1}-\mu_{2}\right)\,f,$
(1)

where $\mu_{i}$ is the chemical potential of the $i$th electronic reservoir
and

$f=\frac{T_{31}-T_{32}}{T_{31}+T_{32}},$ (2)

with $T_{ij}$ the transmission coefficient from Terminal $j$ to Terminal $i$.
Since $f$ gives a measure of the deviation of the voltage drop from the
average value of potentials $\mu_{1}$ and $\mu_{2}$, which is obtained in the
ideal case where the two devices are straight waveguides, we will refer to it
as the _voltage drop deviation_ ; it takes values in the interval [-1, 1].

Figure 1: Schematic representation of the three-terminal system for measuring
the voltage drop through the third terminal located between of two chaotic
devices.

###  II.2 Voltage drop deviation $f$ in terms of scattering elements

Equation (2) is equivalent to the conductance but for three-terminal devices.
Analogously to the conductance the problem of the voltage drop in electronic
devices is reduced to a scattering problem through the quantity $f$. It is
clear that $f$ contains all the information about the system by means of the
transmission coefficients $T_{ij}$. Let us assume that the scattering
properties of each device are known through their scattering matrices, $S_{1}$
for Device 1 and $S_{2}$ for Device 2. The general structure of these matrices
depends on the symmetry properties of the problem, namely

$S_{j}=\left(\begin{array}[]{cc}r_{j}J&t^{\prime}_{j}D_{j}\\\
t_{j}D_{j}^{R}&r^{\prime}_{j}J\end{array}\right),\quad j=1,\,2,$ (3)

where $r_{j}$ ($r^{\prime}_{j}$) and $t_{j}$ ($t^{\prime}_{j}$) are the
reflection and transmission amplitudes for incidence from the left (right) of
the Device $j$. Unitarity is the only requirement on $S_{j}$ in absence of any
symmetry; this case is known as the Unitary symmetry and is labeled by
$\beta=2$ in Dyson’s scheme. In addition, in the presence of time reversal
invariance $S_{j}$ is unitary and symmetric, which corresponds to the
Orthogonal symmetry and is labeled by $\beta=1$; for both cases
$D_{j}=D_{j}^{R}=J=1$. Furthermore, in the presence of time-reversal
invariance but no spin-rotation symmetry, the symmetry of the system is the
Symplectic one, labeled by $\beta=4$, and $S_{j}$ becomes a $4\times 4$ self-
dual matrix RMP69  in which case $J$ is the $2\times 2$ identity matrix,
$J=\openone$, and $D_{j}^{R}$ is the dual matrix of the $2\times 2$ matrix
$D_{j}$, defined by Karol

$D_{j}^{R}=-ZD_{j}^{T}Z,$ (4)

with $D_{j}^{T}$ the transposed matrix of $D_{j}$ and

$Z=\left(\begin{array}[]{cc}0&-1\\\ 1&0\end{array}\right).$ (5)

The dependence of $T_{ij}$ on $S_{1}$ and $S_{2}$ can be obtained explicitly.
The scattering matrix of the three-terminal symmetric configuration is given
by Angel2014

$\mathbf{S}=S_{PP}+S_{PQ}(S_{0}-S_{QQ})^{-1}S_{QP},$ (6)

where $S_{0}$ is the scattering matrix of the junction that accounts for the
coupling to the probe. In fact, $S_{0}$ can be obtained from the experiment
Angel2018 ; Abdu2018 , it reads as

$S_{0}=\frac{1}{3}\left(\begin{array}[]{ccc}-J&2J&2J\\\ 2J&-J&2J\\\
2J&2J&-J\end{array}\right).$ (7)

In Eq. (6) $S_{PP}$ represents the reflections to the terminals, $S_{QP}$ and
$S_{PQ}$ represent the transmissions from the terminals to the inner part of
the system and from the inner part to the terminals; $S_{QQ}$ represents the
internal reflections and $(S_{0}-S_{QQ})^{-1}$ accounts for the multiple
scattering between the Devices and the junction. Explicitly, they are

$\begin{array}[]{cc}\vspace{0.5cm}S_{PP}=\left(\begin{array}[]{ccc}r_{1}J&0&0\\\
0&r^{\prime}_{2}J&0\\\
0&0&0\end{array}\right),&S_{PQ}=\left(\begin{array}[]{ccc}t^{\prime}_{1}D_{1}&0&0\\\
0&t_{2}D^{R}_{2}&0\\\ 0&0&J\end{array}\right),\\\
S_{QP}=\left(\begin{array}[]{ccc}t_{1}D^{R}_{1}&0&0\\\
0&t^{\prime}_{2}D_{2}&0\\\
0&0&J\end{array}\right),&S_{QQ}=\left(\begin{array}[]{ccc}r^{\prime}_{1}J&0&0\\\
0&r_{2}J&0\\\ 0&0&0\end{array}\right).\end{array}$ (8)

By substituting Eqs. (7) and (8) into Eq. (6), the scattering matrix elements
$\mathbf{S}_{31}$ and $\mathbf{S}_{32}$ are obtained, from which
$T_{ij}=|\mathbf{S}_{ij}|^{2}$ for $\beta=1,\,2$, while
$T_{ij}=\frac{1}{2}\mathrm{tr}(\mathbf{S}_{ij}\mathbf{S}^{\dagger}_{ij})$ for
$\beta=4$, and therefore Eq. (2) reduces to

$f=\frac{|t_{1}|^{2}|1+r_{2}|^{2}-|t^{\prime}_{2}|^{2}|1+r^{\prime}_{1}|^{2}}{|t_{1}|^{2}|1+r_{2}|^{2}+|t^{\prime}_{2}|^{2}|1+r^{\prime}_{1}|^{2}}.$
(9)

Note that $S_{31}$, $S_{32}$, and $f$ depend only on the elements of the
individual scattering matrices that describe the devices.

##  III Chaotic scattering for the voltage drop deviation

Of particular interest are the transport properties through chaotic devices,
quantum or classical. The disordered three-terminal system was previously
considered in Ref. Godoy1992b . In any case the scattering quantities, like
the transmission coefficients through each device, fluctuate with respect to a
tuning parameter, like the energy of the incident particles in the quantum
case or the frequency in the classical wave situation, or from sample to
sample. Therefore, it is the distribution of $f$ what becomes more important
rather than a particular value of $f$.

For a chaotic cavity, the statistical fluctuations of the scattering matrix
are described by RMT. There, $S_{j}$ is uniformly distributed according to the
invariant measure $\mathrm{d}\mu_{\beta}(S_{j})$ that defines the circular
ensemble for the symmetry class $\beta$: the Circular Orthogonal Ensemble for
$\beta=1$, the Circular Unitary Ensemble for $\beta=2$, and the Circular
Symplectic Ensemble for $\beta=4$. Hence, the statistical distribution of $f$
can be calculated from the definition

$P_{\beta}(f)=\int\delta\left(f-\frac{T_{31}-T_{32}}{T_{31}+T_{32}}\right)\mathrm{d}\mu_{\beta}(S_{1})\,\mathrm{d}\mu_{\beta}(S_{2}),$
(10)

where $\delta$ is the Dirac delta function.

###  III.1 Statistical distribution of the voltage drop deviation $f$

A useful parameterization for scattering matrices is the polar form

$S_{j}=\left[\begin{array}[]{cc}-\sqrt{1-\tau_{j}}\textnormal{e}^{\textnormal{i}(\phi_{j}+\phi^{\prime}_{j})}J&\sqrt{\tau_{j}}\textnormal{e}^{\textnormal{i}(\phi_{j}+\psi^{\prime}_{j})}D_{j}\\\
\sqrt{\tau_{j}}\textnormal{e}^{\textnormal{i}(\psi_{j}+\phi^{\prime}_{j})}D_{j}^{R}&\sqrt{1-\tau_{j}}\textnormal{e}^{\textnormal{i}(\psi_{j}+\psi^{\prime}_{j})}J\end{array}\right],$
(11)

where $0\leq\tau_{j}\leq 1$ and $\phi_{j}$, $\phi^{\prime}_{j}$, $\psi_{j}$,
$\psi^{\prime}_{j}$ lie in the interval $[0,\,2\pi]$;
$\phi_{j}=\phi^{\prime}_{j}$ and $\psi_{j}=\psi^{\prime}_{j}$ for
$\beta=1,\,4$. Also, for $\beta=4$, $D_{j}$ can be parameterized as

$D_{j}=U_{j}V_{j}^{R}.$ (12)

where $U_{j}$ and $V_{j}$ are $2\times 2$ special unitary matrices and
$V_{j}^{R}$ is the dual of $V_{j}$, as defined in Eq. (4).

For this parameterization, the normalized invariant measure of $S_{j}$ is
given by

$\mathrm{d}\mu_{\beta}(S_{j})=p_{\beta}(\tau_{j})\,\mathrm{d}\tau_{j}\,\frac{\mathrm{d}\phi_{j}}{2\pi}\,\frac{\mathrm{d}\psi_{j}}{2\pi}\times\left\\{\begin{array}[]{lr}\vspace{0.1cm}1\,,&\quad\beta=1,\\\
\vspace{0.1cm}\displaystyle{\frac{\mathrm{d}\phi^{\prime}_{j}}{2\pi}\,\frac{\mathrm{d}\psi^{\prime}_{j}}{2\pi},}&\beta=2\,,\\\
\displaystyle{\mathrm{d}\mu(U_{j})\,\mathrm{d}\mu(V_{j})}\,,&\beta=4,\end{array}\right.$
(13)

where

$p_{\beta}(\tau_{j})=\frac{\beta}{2}\tau_{j}^{\beta/2-1},$ (14)

$\mathrm{d}\mu(U_{j})$ and $\mathrm{d}\mu(V_{j})$ are the invariant measures
of $U_{j}$ and $V_{j}$, respectively. As we will see in what follows, $f$ does
not depend neither on $U_{j}$ nor $V_{j}$ such that we do not need to know
explicitly the expressions for $\mathrm{d}\mu(U_{j})$ and
$\mathrm{d}\mu(V_{j})$ .

In terms of the parameterization of Eq. (11) the voltage drop deviation, given
in Eq. (9), can be written as

$f=1+\frac{\tau_{1}\tau_{2}-2\tau_{2}-2\tau_{2}\sqrt{1-\tau_{1}}\cos{(\psi_{1}+\psi^{\prime}_{1}})}{\tau_{1}+\tau_{2}-\tau_{1}\tau_{2}-\tau_{1}\sqrt{1-\tau_{2}}\cos{(\phi_{2}+\phi^{\prime}_{2})}+\tau_{2}\sqrt{1-\tau_{1}}\cos{(\psi_{1}+\psi^{\prime}_{1}})}.$
(15)

Since $f$ does not depend neither on $U_{j}$ nor $V_{j}$, the integration over
$\mathrm{d}\mu(U_{j})$ and $\mathrm{d}\mu(V_{j})$ in Eq. (10) gives just 1 for
$\beta=4$; except for variables $\tau_{1}$ and $\tau_{2}$, the resulting
expression for $P_{\beta}(f)$ reduces to a one similar to that for $\beta=1$
in the remaining variables. The integration of $\mathrm{d}\phi_{1}/2\pi$,
$\mathrm{d}\phi^{\prime}_{1}/2\pi$, $\mathrm{d}\psi_{2}/2\pi$, and
$\mathrm{d}\psi^{\prime}_{2}/2\pi$ gives 1 for $\beta=2$. Therefore, once we
perform the integration with respect to $\phi_{2}$ and $\phi_{2}^{\prime}$,
making the appropriate change of variables for each symmetry class, and then
with respect to the remaining phases $\psi_{1}$ and $\psi_{1}^{\prime}$, we
arrive at

$P_{\beta}(f)=\int_{0}^{1}\mathrm{d}\tau_{1}\int_{0}^{1}\mathrm{d}\tau_{2}\,p_{\beta}(\tau_{1})\,p_{\beta}(\tau_{2})\,p(f|\tau_{1},\tau_{2}),$
(16)

where $p(f|\tau_{1},\tau_{2})$, which can be interpreted as the conditional
probability distribution of $f$ given $\tau_{1}$ and $\tau_{2}$, is given by

$\displaystyle\ p(f|\tau_{1},\tau_{1})$ $\displaystyle=$
$\displaystyle\frac{1}{\pi^{2}\left(1-f^{2}\right)\sqrt{1-\tau_{1}}}\times$
$\displaystyle\int_{L}\frac{2-\tau_{1}+2\sqrt{1-\tau_{1}}\,x}{\sqrt{\left(\alpha_{+}-x\right)\left(x-\alpha_{-}\right)\left(1-x^{2}\right)}}\,\mathrm{d}x,$

where $L\in(-1,1)\cap(\alpha_{-},\alpha_{+})$ with

$\alpha_{\pm}=\frac{(\tau_{1}-\tau_{2})-(\tau_{1}+\tau_{2}-\tau_{1}\tau_{2})f\pm(1-f)\tau_{1}\sqrt{1-\tau_{2}}}{(1+f)\tau_{2}\sqrt{1-\tau_{1}}},$
(18)

The integration in the domain $L$ is not a trivial problem because it depends
on the variables $\tau_{1}$, $\tau_{2}$, and $f$ through $\alpha_{\pm}$ in a
complicated way, but we can give a further step in the integration by looking
for the limits of integration once the value of $f$ is fixed (see Appendix A
for details). The integral over $x$ in Eq. (III.1) can be expressed in terms
of complete elliptic integrals. The result for $P_{\beta}(f)$ can be expressed
as

$\displaystyle P_{\beta}(f)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi^{2}(1-f^{2})}\times$
$\displaystyle\left[\int_{0}^{1}\mathrm{d}\tau_{1}\frac{p_{\beta}(\tau_{1})}{\sqrt{1-\tau_{1}}}\int_{0}^{u(\tau_{1})}\mathrm{d}\tau_{2}\,p_{\beta}(\tau_{2})\,I_{-1}^{1}(f;\tau_{1},\tau_{2})\right.$
$\displaystyle+$
$\displaystyle\int_{0}^{1}\mathrm{d}\tau_{2}\,p_{\beta}(\tau_{2})\,\int_{0}^{w(\tau_{2})}\mathrm{d}\tau_{1}\frac{p_{\beta}(\tau_{1})}{\sqrt{1-\tau_{1}}}\,I_{\alpha_{-}}^{\alpha_{+}}(f;\tau_{1},\tau_{2})$
$\displaystyle+$ $\displaystyle
H(f)\,H(1-f)\int_{0}^{1}\mathrm{d}\tau_{1}\frac{p_{\beta}(\tau_{1})}{\sqrt{1-\tau_{1}}}\times$
$\displaystyle\int_{u(\tau_{1})}^{v(\tau_{1})}\,\mathrm{d}\tau_{2}\,p_{\beta}(\tau_{2})\,I_{-1}^{\alpha_{+}}(f;\tau_{1},\tau_{2})$
$\displaystyle+$
$\displaystyle\left.H(-f)\,H(f+1)\int_{0}^{1}\mathrm{d}\tau_{1}\frac{p_{\beta}(\tau_{1})}{\sqrt{1-\tau_{1}}}\right.\times$
$\displaystyle\left.\int_{u(\tau_{1})}^{v(\tau_{1})}\mathrm{d}\tau_{2}\,p_{\beta}(\tau_{2})\,I_{\alpha_{-}}^{1}(f;\tau_{1},\tau_{2})\right],$

where $H(x)$ is the Heaviside step function and the limits are

$\displaystyle u(\tau_{1})$ $\displaystyle=$
$\displaystyle\frac{r(2s_{-}-r)}{s_{-}^{2}},\quad
v(\tau_{1})=\frac{r(2s_{+}-r)}{s_{+}^{2}},\quad$ (20)
$\displaystyle\mathrm{and}\quad
w(\tau_{2})=\frac{r^{\prime}(2s^{\prime}-r^{\prime})}{{s^{\prime}}^{2}},$

where $r^{\prime}=r(\tau_{1}\to\tau_{2})$ and
$s^{\prime}=s_{-}(\tau_{1}\to\tau_{2})$ with

$\displaystyle r(\tau_{1})$ $\displaystyle=$
$\displaystyle(1-f)\tau_{1}\quad\mbox{and}\quad$ $\displaystyle
s_{\pm}(\tau_{1})$ $\displaystyle=$ $\displaystyle
1+(1-\tau_{1})f\pm(1+f)\sqrt{1-\tau_{1}}.$ (21)

In Eq. (III.1),

$\displaystyle
I_{-1}^{1}(f;\tau_{1},\tau_{2})=\frac{2}{\sqrt{(\alpha_{+}+1)(1-\alpha_{-})}}$
$\displaystyle\times$
$\displaystyle[(A+B\alpha_{-})K(k_{1})-(1+\alpha_{-})B\,\Pi(r_{1}^{2},k_{1})],$
$\displaystyle
I_{\alpha_{-}}^{\alpha_{+}}(f;\tau_{1},\tau_{2})=\frac{2}{\sqrt{(\alpha_{+}+1)(1-\alpha_{-})}}$
$\displaystyle\times$
$\displaystyle[(A-B)K(k_{1})+(1+\alpha_{-})B\,\Pi({r^{\prime}_{1}}^{2},k_{1})],$
$\displaystyle
I_{-1}^{\alpha_{+}}(f;\tau_{1},\tau_{2})=\frac{\sqrt{2}}{\sqrt{\alpha_{+}-\alpha_{-}}}$
$\displaystyle\times$
$\displaystyle[(A+B\alpha_{-})K(k_{2})-(1+\alpha_{-})B\,\Pi({r^{\prime}_{2}}^{2},k_{2})],$

and

$\displaystyle
I_{\alpha_{-}}^{1}(f;\tau_{1},\tau_{2})=\frac{\sqrt{2}}{\sqrt{\alpha_{+}-\alpha_{-}}}$
$\displaystyle\times$
$\displaystyle[(A-B)K(k_{2})+(1+\alpha_{-})B\,\Pi({r_{2}}^{2},k_{2})],$

where $A=2-\tau_{1}$ and $B=2\sqrt{1-\tau_{1}}$; $K(k)$ and $\Pi(r^{2},k)$ are
complete elliptic integrals of the first and third kind, respectively, and

$\displaystyle k_{1}^{2}$ $\displaystyle=$
$\displaystyle\frac{2(\alpha_{+}-\alpha_{-})}{(\alpha_{+}+1)(1-\alpha_{-})},\quad
r_{1}^{2}=\frac{2}{1-\alpha_{-}},$ (26)
$\displaystyle\quad\mbox{and}\quad{r^{\prime}_{1}}^{2}=\frac{\alpha_{+}-\alpha_{-}}{\alpha_{+}+1};$
$k_{2}^{2}=\frac{1}{k_{1}^{2}},\quad
r_{2}^{2}=\frac{1}{r_{1}^{2}},\quad\mbox{and}\quad{r^{\prime}_{2}}^{2}=\frac{1}{{r^{\prime}_{1}}^{2}}.$
(27)

Figure 2: Comparison between the analytical results (continuous lines) with
random matrix simulations (histograms) for $P_{\beta}(f)$: (a) $\beta=1$, (b)
$\beta=2$, and (c) $\beta=4$.

The remaining integrals with respect to $\tau_{1}$ and $\tau_{2}$ can be
performed numerically. In order to verify our results, in Fig. 2 we compare
them with the statistical distributions obtained from numerical simulations
using Eqs. (11)-(15). An excellent agreement is observed.

We can observe in Fig. 2 that the distribution of $f$ is symmetric with
respect to $f=0$. Moreover, it diverges at $f=\pm 1$ for $\beta=1$, while it
shows finite peaks at $f=\pm 1$ for $\beta=2$, which is a clear effect of the
weak localization. For $\beta=4$ the distribution becomes zero at $f=\pm 1$
due to the antilocalization phenomenon. The differences found in the
$P_{\beta}(f)$ between the symmetry classes are important signatures of the
chaotic setup we consider here, since in the equivalent disordered
configuration no differences were found Godoy1992b .

##  IV Experimental realizations with microwave graphs

Figure 3: (a) Photograph of the three-port experimental setup for the case of
time-reversal symmetry ($\beta=1$). For the realization of a break of time-
reversal symmetry ($\beta=2$) one T junction in each of the subgraphs is
replaced by a circulator. Both subgraphs are geometrically different but have
the same total length, i. e., the mean level density is the same. (b) Sketch
of the graph for $\beta=1$ (GOE) and $\beta=2$ (GUE), where for $\beta=1$ the
circulators are replaced by ordinary T junctions (red). For a realization of
$\beta=4$ (GSE) each GUE subgraph is complemented by another geometrically
identical one but with an opposite sense of rotation of the circulators
(grey). The two respective copies are connected by pairs of bonds with length
differences corresponding to a phase difference of $\pi$ for the propagating
waves (see Ref. Abdu2016  for details, where also a photograph of a single GSE
graph can be found).

Since the voltage drop deviation depends only on the scattering properties of
the devices, it can be emulated in classical wave systems, in particular
experiments with microwave graphs can be performed. Figure 3(a) shows a
photograph of the experimental set-up for the case of time-reversal invariance
($\beta=1$). The devices consist of chaotic microwave networks formed by
coaxial semirigid cables (Huber & Suhner EZ-141) with SMA connectors, coupled
by T junctions at the nodes. One microwave port attached to the left subgraph
acts as the input, another one on the right as the output, and a third port
attached to the connecting cable between the two subgraphs as the probe. To
realize a break of time-reversal invariance ($\beta=2$), in each of the two
subgraphs one of the T junctions is replaced by a circulator (Aerotek
I70-1FFF) with an operating frequency range from 6 to 12 GHz. A circulator
introduces directionality, waves entering via port 1 leave via port 2. The
transmission intensities $T_{31}$ and $T_{32}$ were measured by an Agilent
8720ES vector network analyzer (VNA). The $\beta=1$ case has been additionally
realized in a billiard setup. Details are presented in Appendix B.

For the realization of the $\beta=4$ (GSE) case two GSE graphs are needed.
Since each GSE graph is composed by two GUE subgraphs, representing the spin-
up and spin-down components Abdu2016 , we need a total of four subgraphs.
Figure 3(b) shows a sketch of the three-terminal GSE graph used to measure the
scattering matrix. An essential ingredient of the setup are two pairs of bonds
with length differences $\Delta l$ corresponding to phase differences of
$\Delta\varphi=\pi$ for the propagating waves. In the experiment we took
spectra for fixed $\Delta l$ and converted them into spectra for fixed
$\Delta\varphi$ using $\Delta\varphi=k\Delta l$, where $k$ is the wavenumber.
Details can be found in Ref. Abdu2016 . The ports now appear in pairs
$1,\bar{1},2,\bar{2},3,\bar{3}$, and the scattering matrix elements turn into
$2\times 2$ matrices

$\displaystyle\mathbf{S}_{ij}=\left(\begin{array}[]{cc}S_{ij}&S_{i\bar{j}}\\\
S_{\bar{i}j}&S_{\bar{i}\bar{j}}\end{array}\right).$ (30)

In the spirit of the spin analogy, $S_{ij}$ and $S_{i\bar{j}}$ correspond to
transmissions without and with spin-flip, respectively. $\mathbf{S}_{ij}$ may
be written in terms of quaternions as
$\mathbf{S}_{ij}=S^{0}_{ij}\openone+\mathrm{i}S^{x}_{ij}\sigma_{x}+\mathrm{i}S^{y}_{ij}\sigma_{y}+\mathrm{i}S^{z}_{ij}\sigma_{z}$,
where for a symplectic symmetry all coefficients $S^{0}_{ij}$, $S^{x}_{ij}$,
etc. are real numbers and $\sigma_{j}$ is the corresponding Pali matrix $j$.
As a consequence $\mathbf{S}_{ij}\mathbf{S}_{ij}^{\dagger}$ is a multiple of
the unit matrix which allows for a simple check of the quality of the
realization of the setup with symplectic symmetry. In our experiments
$\mathbf{S}_{31}$ and $\mathbf{S}_{32}$ were found to be quaternion real
within an error of $2.19\%$ and $1.93\%$, respectively. The mentioned multiple
is nothing but the transmission coefficient, hence
$T_{31}=\frac{1}{2}\mathrm{tr}(\mathbf{S}_{31}\mathbf{S}_{31}^{\dagger})$ and
$T_{32}=\frac{1}{2}\mathrm{tr}(\mathbf{S}_{32}\mathbf{S}_{32}^{\dagger})$.

Figure 4: (Color online) Experimental voltage drop deviation $f$ as a function
of frequency is shown in the upper panels for the three symmetry classes.
Their corresponding statistical distribution are shown as histograms. The
continuous (blue) lines in the lower panels correspond to the theoretical
result shown in Fig 2, while the dashed (orange) lines are the RMT simulations
which take into account the dissipation and imperfect coupling: $\beta=1$
(left), $\beta=2$ (middle), and $\beta=4$ (right). For comparison purposes, in
the insets in the lower panels we show the difference between the numerical
and the experimental distribution $\delta
P_{\beta}(f)=P_{\beta}(f)_{\mathrm{num}}-P_{\beta}(f)_{\mathrm{expt}}$. For
the statistical analysis we used an ensemble of $5\times 10^{4}$ realizations.

The experimental distribution of $f$ is shown in Fig. 4 as histograms for the
three symmetry classes $\beta$, where it is compared with the theoretical
result given by Eq. (III.1). As can be observed there is a qualitatively good
agreement for all $\beta$. Despite the fact that there is a quantitative
difference between theory and experiment, which is due to the phenomena of
dissipation and imperfect coupling between the graphs and ports, this
difference is not so large. In Fig. 4 we also show the corrected distribution
obtained from RMT simulations, using the Heidelberg approach (see below), once
these two phenomena, dissipation and imperfect coupling, are taken into
account.

The dissipation and imperfect coupling can be quantified by two parameters:
$T_{1}$ for the coupling strength and $\gamma$ for the dissipation. The
coupling between the graphs and ports $T_{1}$ is extracted from the
experimental data as $T_{1}=1-|\langle S_{11}\rangle|^{2}$, where $\langle
S_{11}\rangle$ is the average (with respect to the frequency) of the
scattering matrix element 11. Here we obtained $T_{1}=0.99$ for $\beta=1$,
$T_{1}=0.99$ for $\beta=2$, and $T_{1}=0.953$ for $\beta=4$. As can be seen
the coupling strength is almost 1, which means that the coupling between
graphs and ports, although not perfect, can be considered as very good. The
effect of dissipation is clearly observed in the experimental transmission
intensities $T_{31}$ and $T_{32}$. Their fluctuations as a function of the
frequency are shown in Fig. 5, where we observe that they do not reach the
value 1. The dissipation parameter $\gamma$ can be quantified by fitting the
autocorrelation function $C_{11}(t)$ of the 11 element of the scattering
matrix Schafer2003 ; Angel2018 ; Fyodorov2005b

$\frac{C_{11}(t)}{T_{1}^{2}}=\left\\{\begin{array}[]{ll}\vspace{0.1cm}\left[\displaystyle{\frac{3}{(1+2T_{1}t)^{3}}-\frac{b_{1,2}(t)}{(1+T_{1}t)^{4}}}\right]\mathrm{e}^{-\gamma
t}&\mbox{for}\,\,\beta=1,\\\
\vspace{0.1cm}\left[\displaystyle{\frac{2}{(1+T_{1}t)^{4}}-\frac{2^{6}b_{2,2}(t)}{(2+T_{1}t)^{6}}}\right]\mathrm{e}^{-\gamma
t}&\mbox{for}\,\,\beta=2,\\\
\left[\displaystyle{\frac{6}{(1+T_{1}t)^{6}}-\frac{2^{12}b_{4,2}(t)}{(2+T_{1}t)^{10}}}\right]\mathrm{e}^{-2\gamma
t}&\mbox{for}\,\,\beta=4,\end{array}\right.$ (31)

where $b_{\beta,2}(t)$ is the two-level form factor Guhr1998  and $T_{1}$,
$\gamma$, and $t$ are given in dimensionless units. The best fit achieved for
$C_{11}(t)$, for $t<1$, is shown in Fig. 5 for all symmetry classes, for which
we obtain $\gamma=2.6$ for $\beta=1$, $\gamma=2.0$ for $\beta=2$, and
$\gamma=1.8$ for $\beta=4$.

Figure 5: Experimental transmission intensities $T_{31}$ and $T_{32}$ as a
function of frequency are shown in upper and middle panels for the three
symmetry classes. The lower panels show the autocorrelation function: the
fluctuations are the experimental measurements and the dashed lines are the
best fits of Eq. (31) to the data for $t<1$.

Once these parameters are determined, RMT simulations can be performed for the
scattering matrix of each graph, assumed as quantum systems. In the Heidelberg
approach the scattering matrix of a graph can be generated as Brouwer1997

$S_{g}=1-2\pi\mathrm{i}W^{\dagger}\frac{1}{E-\mathcal{H}+\mathrm{i}\pi
WW^{\dagger}}W,$ (32)

where $E$ represents the energy of the incoming wave and $W$ is the matrix
that couples the open modes in the ports to the internal modes of the system.
$\mathcal{H}$ is an effective Hamiltonian that includes the dissipation
$\gamma$,
$\mathcal{H_{\mu\nu}}=H_{\mu\nu}-\mathrm{i}\delta_{\mu\nu}\gamma\Delta/4\pi$,
with $\Delta$ the mean level spacing. The imperfect coupling can be modeled by
adding identical barriers, with transmission intensity $T_{1}$, between the
graph 1 and port 1, between the graphs and the T junction, and between graph 2
and port 2. The numerical simulations of the three-terminal system, with the
obtained values of $T_{1}$ and $\gamma$, lead us to the result shown in Fig. 4
for the distribution of $f$; see the orange dashed lines in the lower panels.
As can be observed, the agreement is very good.

##  V Conclusions

A three-terminal symmetrical system consisting of two microwave graphs was
used to measure indirectly the voltage drop between the graphs. One port was
used as an input, a second port as an exit, and a third port as a probe. The
fluctuations of the quantity $f$, that accounts for the deviation from the
mean value of the potentials, are explained qualitatively by an ideal result
obtained from the scattering approach of random matrix theory for the three
Dyson’s symmetry classes. We found that the distribution of $f$ is symmetric
with respect to zero, whose shape is quite similar to the one obtained in the
disordered case, in the insulating regime, but with an important difference
which is the effect of weak localization and antilocalization phenomena not
found in the disordered case. A more accurate description is obtained when
dissipation and imperfect coupling between the ports and graphs are taken into
account.

###### Acknowledgements.

This work was supported by CONACyT (Grant No. CB-2016/285776). F. C.-R. thanks
financial support from CONACyT and A.M.M-A. acknowledges support from DGAPA-
UNAM. J.A.M.-B. acknowledges financial support from CONACyT (Grant No.
A1-S-22706). The experiments were funded by the Deutsche
Forschungsgemeinschaft via the individual grants STO 157/17-1 and KU 1525/3-1
including a short-term visit of A.M.M.-A. in Marburg.

##  Appendix A Limits of integration in Eq. (III.1)

The integration with respect to $x\in(-1,1)$ in Eq. (III.1) should be
performed in the interval $L$, which is the intersection of the intervals
$(-1,1)$ and $(\alpha_{-},\alpha_{+})$. Since $\alpha_{+}$ and $\alpha_{-}$
depend on the values of $\tau_{1}$, $\tau_{2}$, and $f$, the values of
$\tau_{1}$ and $\tau_{2}$ may be affected by that intersection for fixed $f$.
This leads to integrals of the form

$I_{a}^{b}(f,\tau_{1},\tau_{2})=\int_{a}^{b}\frac{(A+Bx)dx}{\sqrt{(\alpha_{+}-x)(x-\alpha_{-})(1+x)(1-x)}},$
(33)

which gives rise to complete elliptic integrals of the first and third kind.

Figure 6: (a) Intervals of integration for $x$ which imply restrictions on
$\tau_{1}$ and $\tau_{2}$ (shaded regions): for $|f|=0.6$ (b)
$\tau_{2}\in(0,u(\tau_{1}))$ with $\tau_{1}\in(0,1)$; (c)
$\tau_{1}\in(0,w(\tau_{2}))$ with $\tau_{2}\in(0,1)$; (d)
$\tau_{2}\in(u(\tau_{1}),v(\tau_{1}))$ for $f>0$ and
$\tau_{2}\in(v(\tau_{1}),u(\tau_{1}))$ for $f<0$, with $\tau_{1}\in(0,1)$. The
functions $u(\tau)$, $v(\tau)$, and $w(\tau)$ are defined in Eq. (20).

Four conditions arise, as illustrated in Fig. 6(a). Restrictions on $\tau_{1}$
and $\tau_{2}$ are obtained for a fixed value of $f$. These restrictions lead
to the several regions in the plane $\tau_{1}\tau_{2}$, as can be seen in Fig.
6.

#### Condition $C_{1}$

$\alpha_{-}\leq-1$ and $1\leq\alpha_{+}$. Under these conditions $x$ runs over
its full domain, $x\in(-1,1)$. Restriction on $\tau_{2}$ for a fixed value of
$f$ is obtained from these two conditions. The first condition,
$\alpha_{-}\leq-1$, leads to $\tau_{2}\geq 0$ while the second condition,
$1\leq\alpha_{+}$, restricts $\tau_{2}$ to be smaller than $u(\tau_{1})$,
where $u(\tau_{1})$ is given in equation of (20). Therefore,
$\tau_{2}\in(0,u(\tau_{1}))$ with $\tau_{1}\in(0,1)$.

#### Condition $C_{2}$

$-1\leq\alpha_{-}\leq\alpha_{+}\leq 1$. For this condition
$x\in(\alpha_{-},\alpha_{+})$. Here, it is more convenient to see the
restriction on $\tau_{1}$ for which we obtain that
$\tau_{1}\in(0,w(\tau_{2}))$, where $w(\tau_{2})$ is given by the third
equation in (20), for $\tau_{2}\in(0,1)$.

#### Condition $C_{3}$

$\alpha_{-}\leq-1$ and $-1\leq\alpha_{+}\leq 1$. The integration over $x$ is
defined in the interval $(-1,\alpha_{+})$. The first part of the condition is
the same as that for Condition $C_{1}$, which restricts to $\tau_{2}\geq 0$.
The second condition has two parts, a first one is $\alpha_{+}\leq 1$ whose
result is similar to the corresponding second part in Condition $C_{1}$, but
the inequality in opposite sense; hence, $\tau_{2}\geq u(\tau_{1})$. In a
similar manner, the second part $\alpha_{+}\geq-1$ leads to $\tau_{2}\leq
v(\tau_{1})$ for $f>0$, which is given by the second equation in (20).
Therefore, $\tau_{2}\in(u(\tau_{1}),v(\tau_{1}))$ with $\tau_{1}\in(0,1)$ and
$f>0$.

#### Condition $C_{4}$

$-1\leq\alpha_{-}\leq 1$ and $1\leq\alpha_{+}$. The integration over $x$ is
defined in the interval $(\alpha_{-},1)$. The second part of the condition is
the same as that for Condition $C_{1}$ which restricts to
$\tau_{2}\in(0,u(\tau_{1}))$ if $f$ is positive. The first condition has two
parts, the first one, $\alpha_{-}\leq 1$ leads to the same restriction to
$\tau_{2}$. The second part gives $\tau_{2}\in(v(\tau_{1}),1)$ for $f>0$. This
result is equivalent to $\tau_{2}\in(u(\tau_{1}),v(\tau_{1}))$ for $f<0$.

To apply these conditions to the integral of Eq. (33) we use the equations
254.00, 254.10, 336.01, 336.60, and 340.04 of Ref. ByrdHandbook , then we
integrate with respect to $\tau_{1}$ and $\tau_{2}$, and finally we arrive at
Eq. (III.1).

##  Appendix B Three-port microwave billiard experiment

Figure 7: (a) Sketch of the three-port microwave billiard in the presence of
$\beta=1$ symmetry. The ports are labeled as 1, 2, and 3, where port 3 is used
as a probe. The experimental transmission intensities $T_{31}$ and $T_{32}$
are measured between ports 1 and 3, and between ports 2 and 3, respectively.
(b) Distribution $P_{1}(f)$: the continuous (black) line corresponds to the
theoretical result shown in Fig. 2, the dashed (red) line is the RMT
simulation which takes into account dissipation and imperfect coupling, while
the histogram corresponds to the experimental result.

For the $\beta=1$ case, in addition to the microwave graph experiments an
experiment in a billiard setup has been performed. A sketch is shown in Fig.
7(a). The billiard is constructed on an aluminum plate with two sub-billiards
of the same shape separated by a central bar. The mirror symmetry is broken by
two semicircular obstacles, attached to the bottom boundary. The ports are
labeled as 1, 2, and 3, where port 3 is used as a probe. The experimental
transmission intensities $T_{31}$ and $T_{32}$ have been measured between port
1 and port 3, and between port 2 and port 3, respectively, for frequencies
from 1 to 17 GHz. With a distance of $d=8$ mm between top and bottom plate the
billiard is quasi-two-dimensional in the whole frequency range.

In panel (b) of Fig. 7 we show as histograms the experimental distribution of
$P_{1}(f)$ obtained from the graph (black) and the billiard (blue) setting.
The analytical result is shown in the continuous (black) line, while the RMT
simulation, which take into account the dissipation and imperfect coupling, is
shown in (red) dashed line. Again, a good agreement between experiment and
theory is found for $P_{1}(f)$.

## References

  * (1) Y. V. Fyodorov, T. Kottos, H.-J. Stöckmann, Trends in quantum chaotic scattering, J. Phys. A: Math. Gen. 38 (2005) Preface.
  * (2) H. Schanze, E. R. P. Alves, C. H. Lewenkopf, and H.-J. Stöckmann, Phys. Rev. E 64, 065201(R) (2001).
  * (3) H. Schanze, H.-J. Stöckmann, M. Martínez-Mares, and C. H. Lewenkopf, Phys. Rev. E 71, 016223 (2005).
  * (4) R. Schäfer, T. Gorin, T. H. Seligman, and H.-J. Stöckmann, J. Phys. A: Math. Gen. 36, 3289 (2003).
  * (5) U. Kuhl, M. Martínez-Mares, R. A. Méndez-Sánchez, and H.-J. Stöckmann, Phys. Rev. Lett. 94, 144101 (2005).
  * (6) M. Büttiker, IBM J. Res. Develop. 32, 317 (1988).
  * (7) K. Schwab, E. A. H. J. M. Worlock, and M. L. Roukes, Nature 404, 974 (2000).
  * (8) M. Büttiker, Y. Imry, and M. Y. Azbel, Phys. Rev. A 30, 1982 (1984).
  * (9) P. W. Brouwer and C. W. J. Beenakker, Phys. Rev. B 55, 4695 (1997).
  * (10) M. W. Keller, A. Mittal, J. W. Sleight, R. G. Wheeler, D. E. Prober, R. N. Sacks, and H. Shtrikmann, Phys. Rev. B 53, R1693 (1996).
  * (11) C. M. Marcus, A. J. Rimberg, R. M. Westervelt, P. F. Hopkins, and A. C. Gossard, Phys. Rev. Lett. 69, 506 (1992) .
  * (12) I. H. Chan, R. M. Clarke, C. M. Marcus, K. Campman, and A. C. Gossard, Phys. Rev. Lett. 74, 3876 (1995).
  * (13) M. Martínez-Mares, Phys. Rev. E 72, 036202 (2005).
  * (14) E. Flores-Olmedo, A. M. Martínez-Argüello, M. Martínez-Mares, G. Báez, J. A. Franco-Villafañe, and R. A. Méndez-Sánchez, Sci. Rep. 6, 25157 (2016).
  * (15) M. Büttiker, Four-terminal phase-coherent conductance, Phys. Rev. Lett. 57, 1761 (1986).
  * (16) V. A. Gopar, M. Martínez, and P. A. Mello, Phys. Rev. B 50, 2502 (1994).
  * (17) L. Arrachea, C. Naón, and M. Salvay, Phys. Rev. B 77, 233105 (2008).
  * (18) C. Texier and G. Montambaux, Physica E 82, 272 (2016).
  * (19) F. Foieri, L. Arrachea, and M. J. Sánchez, Phys. Rev. B 79, 085430 (2009).
  * (20) J. L. D’Amato and H. M. Pastawski, Phys. Rev. B 41, 7411 (1990).
  * (21) C. J. Cattena, L. J. Fernández-Alcázar, R. A. Bustos-Marún, D. Nozaki, and H. M. Pastawski, J. Phys.: Condens. Matter 26, 345304 (2014).
  * (22) A. M. Song, A. Lorke, A. Kriele, J. P. Kotthaus, W. Wegscheider, and M. Bichler, Phys. Rev. Lett. 80, 3831 (1998).
  * (23) B. Gao, Y. F. Chen, M. S. Fuhrer, D. C. Glattli, and A. Bachtold, Phys. Rev. Lett. 95, 196802 (2005).
  * (24) A. M. Martínez-Argüello, A. Rehemanjiang, M. Martínez-Mares, J. A. Méndez-Bermúdez, H.-J. Stöckmann, and U. Kuhl, Phys. Rev. B 98, 075311 (2018).
  * (25) A. M. Martínez-Argüello, J. A. Méndez-Bermúdez, and M. Martínez-Mares, Phys. Rev. E 99, 062202 (2019).
  * (26) S. Godoy and P. A. Mello, EPL 17, 243 (1992).
  * (27) S. Godoy and P. A. Mello, Phys. Rev. B 46, 2346 (1992).
  * (28) A. M. Martínez-Argüello, E. Castaño, and M. Martínez-Mares, Random matrix study for a three-terminal chaotic device, in _Special Topics on Transport Theory: Electrons, Waves, and Diffusion in Confined Systems_ , AIP Conf. Proc. No. 1579 (AIP, Melville, NY , 2014), p. 46.
  * (29) C. W. J. Beenakker, Rev. Mod. Phys. 69, 731 (1997).
  * (30) K. Życzkowski, _Random Matrices of Circular Symplectic Ensemble_ in _Chaos–The Interplay Between Stochastic and Deterministic Behaviour_ , Lecture Notes in Physics, vol 457, edited by P. Garbaczewski, M. Wolf, and A. Weron (Springer, Berlin, Heidelberg, 1995).
  * (31) A. Rehemanjiang, M. Richter, U. Kuhl, and H.-J. Stöckmann, Phys. Rev. E 97, 022204 (2018).
  * (32) A. Rehemanjiang, M. Allgaier, C. H. Joyner, S. Müller, M. Sieber, U. Kuhl, and H.-J. Stöckmann, Phys. Rev. Lett. 117, 064101 (2016).
  * (33) Y. V. Fyodorov, D. V. Savin, and H.-J. Sommers, J. Phys. A: Math. Gen. 38, 10731 (2005).
  * (34) T. Guhr, A. Müller-Groeling, and H. A. Weidenmüller, Phys. Rep. 299, 189 (1998).
  * (35) P. F. Byrd and M. D. Friedman, Handbook of Elliptic Integrals for Engineers and Physicists, Springer-Verlag Berlin Heidelberg GMBH, Berlin, Heidelberg, 1954.


# DHR: Distributed Hybrid Rendering for Metaverse Experiences

Tan Yu Wei  National University of SingaporeSingapore<EMAIL_ADDRESS>,  Alden
Tan  National University of SingaporeSingapore<EMAIL_ADDRESS>,  Nicholas
Nge  National University of SingaporeSingapore<EMAIL_ADDRESS>and  Anand
Bhojan  National University of SingaporeSingapore<EMAIL_ADDRESS>
(2022)

###### Abstract.

Classically, rasterization techniques are performed for real-time rendering to
meet the constraint of interactive frame rates. However, such techniques do
not produce realistic results as compared to ray tracing approaches. Hence,
hybrid rendering has emerged to improve the graphics fidelity of rasterization
with ray tracing in real-time. We explore the approach of distributed
rendering in incorporating real-time hybrid rendering into metaverse
experiences for immersive graphics. In standalone extended reality (XR)
devices, such ray tracing-enabled graphics is only feasible through pure
cloud-based remote rendering systems that rely on low-latency networks to
transmit real-time ray-traced data in response to interactive user input.
Under high network latency conditions, remote rendering might not be able to
maintain interactive frame rates for the client, adversely affecting the user
experience. We adopt hybrid rendering via a distributed rendering approach by
integrating ray tracing on powerful remote hardware with raster-based
rendering on user access devices. With this hybrid approach, our technique can
help standalone XR devices achieve ray tracing-incorporated graphics and
maintain interactive frame rates even under high-latency conditions.

real-time, ray tracing, hybrid rendering, distributed rendering, interactive
applications, metaverse

††journalyear: 2022††copyright: rightsretained††conference: Proceedings of the
1st Workshop on Interactive eXtended Reality; October 14, 2022; Lisboa,
Portugal††booktitle: Proceedings of the 1st Workshop on Interactive eXtended
Reality (IXR ’22), Oct. 14, 2022, Lisboa, Portugal††doi:
10.1145/3552483.3556455††isbn: 978-1-4503-9501-4/22/10††submissionid:
ixr008††ccs: Computing methodologies Rendering††ccs: Computing methodologies
Ray tracing

##  1\. Introduction

The metaverse refers to a collective virtual environment on the Internet that
incorporates both the physical and digital world (Lee et al., 2021) as
empowered by extended reality (XR) technology for immersive real-time
interaction and collaboration. According to virtual reality (VR) pioneer and
metaverse expert Tony Parisi (Parisi, 2021), access devices for the metaverse
are not limited to standalone XR devices but can also include mobile phones
and display walls. However, he acknowledges that one’s experience of the
metaverse can be enhanced with “immersive hardware devices”. Hence, we target
standalone XR devices to elevate the quality of the metaverse.

Ray tracing is a rendering technique which generates realistic lighting and
camera effects by simulating the propagation of light with respect to scene
geometry and the camera. However, it is computationally heavy due to the costs
of calculating ray-geometry intersections and multiple shading per pixel. In
contrast, rasterization provides a fast approximation of ray tracing’s visual
quality but suffers from certain inaccuracies which cannot be solved with just
screen-space information. Given the advancements in hardware-accelerated
personal computer (PC) graphics processing units (GPUs) that support ray
tracing such as the NVIDIA GeForce RTX, we can now achieve more convincing
graphics in real-time with a reasonable additional performance cost using ray
tracing techniques that were once only feasible in offline rendering.

However, for standalone XR devices such as VR headsets, GPUs that support
hardware-accelerated ray tracing are not feasible to realise with current
technologies due to their size factor, heat generation and power consumption.
Hence, such system hardware cannot perform real-time ray tracing fast enough
to meet the requirements of 90 frames per second (FPS) for interactive XR
applications. Nonetheless, ray tracing-based rendering can generate advanced
lighting and camera effects for more photorealistic graphics as compared to
traditional rasterization, helping to simulate a more immersive digital
environment for the metaverse. Hybrid rendering for real-time applications
seeks to combine ray tracing and rasterization techniques for a performance-
quality tradeoff, producing higher fidelity results than classic real-time
rendering while still maintaining interactive frame rates. However, on
standalone XR devices, even hybrid rendering with only some ray tracing is not
feasible to meet the tight performance constraint.

In cloud gaming, remote rendering is often performed where the user’s access
device (the client) only handles the player’s inputs and sends them to the
remote server for rendering. After the rendering is done, output video frames
are then streamed to the client to be displayed to the user. Hence, hybrid
rendering can be performed as part of remote rendering which leverages the
graphics capability of remote servers for rendering as compared to local
rendering within the device itself. However, the cloud gaming industry is
still facing multiple challenges including well-known latency and bandwidth
issues. For instance, latency limits the proximity that a cloud server can
support and in turn the number of clients, making cloud gaming not as
attractive as other cloud-based services. Hence, the industry is focusing on
cloud-assisted rather than cloud-dependent systems.

We propose the novel cloud-assisted approach of hybrid rendering distributed
between the cloud system and the client which we call Distributed Hybrid
Rendering (DHR). DHR leverages both the client device and remote server for
rendering by allocating the rendering workload based on real-time ray tracing
capability. The client and server work together to generate ray tracing-
incorporated high-fidelity graphics while maintaining interactive frame rates
for metaverse experiences. By distributing the rendering workload, we can
leverage the client devices’ limited capability for rasterization to also
serve as a backup should the ray-traced data from the remote server be delayed
or lost. When this data is missing, the use of the application is not
disrupted as local rasterization in the client with some approximation
mechanisms can still provide desirable visual quality. DHR can be added
directly to the rendering engine with no interference from the application
developer.

This paper is an extended version of our previous work Tan et al. (2021) with
the following key contributions:

  * •

Details of the design and implementation of a basic DHR prototype for dynamic
ray-traced shadows.

  * •

A frame prediction mechanism to maintain spatial detail and temporal coherence
under adverse network conditions.

  * •

Qualitative and quantitative evaluation of our technique, including the
analysis of key factors such as network ping and the number of predicted
frames.

##  2\. Related Work

Hybrid rendering combines ray tracing and rasterization to generate high-
quality graphics while maintaining interactive frame rates for real-time
applications. Instead of full ray tracing which might be too computationally
expensive for geometrically complex scenes, ray tracing is typically at most
partially employed such as for certain lighting effects (Cabeleira, 2010;
Barré-Brisebois et al., 2019) or pixels (Marrs et al., 2018; Macedo et al.,
2018; Beck et al., 1981; Hertel et al., 2009; Lauterbach and Manocha, 2009)
where the visual result of rasterization is less desirable. We employ hybrid
rendering to produce ray tracing-enabled high-fidelity graphics for metaverse
experiences while meeting interactive performance budgets.

Among the four main delivery modes for game contents (Anand and Wenren, 2017),
video streaming and graphics streaming (or image-based streaming and
instruction-based streaming as termed by Chan et al. (2017)) are typically
used with regards to rendering. Image-based streaming is where the server
performs the rendering (Holthe et al., 2009), and instruction-based streaming
is where the server only computes the graphics commands for the rendering that
the client proceeds to perform (Eisert and Fechteler, 2008). Although
instruction-based streaming reduces the amount of data that needs to be
transmitted (i.e. graphics commands instead of whole video frames), image-
based streaming is the convention for cloud gaming as rendering can be
performed on powerful remote GPUs to achieve high-quality graphics even for
thin-client user access devices. Popular cloud gaming services like Google
Stadia (Google, 2021), NVIDIA GeForce NOW (NVIDIA, 2021) and Amazon Luna
(Amazon, 2021), all employ image-based streaming (i.e. remote rendering).
Standalone XR devices are able to produce realistic visual results through the
remote rendering of hybrid rendering techniques, provided that the network
latency is low for the interactive streaming of video content to the client.
However, under high-latency conditions, this pure cloud-based solution may not
fulfil the requirement of interactive frame rates (Anand and Hao Edwin, 2014;
Anand and Wenren, 2017) if data transmission is too slow, resulting in
undesired lag for the user.

To reap the benefits of both the image- and instruction-based streaming modes,
Chan et al. (2017) propose a hybrid-streaming workflow that uses image-based
streaming for far objects and instruction-based streaming for near objects,
overlaying the results in order of depth. This workflow adopts the approach of
collaborative or distributed rendering (Cuervo et al., 2015) where rendering
is performed on both the client and server and their results are combined for
the final output. Other distributed rendering techniques include layered
coding (Chuah and Cheung, 2014; Chuah et al., 2016) where the server renders
two versions of the frame of different quality based on polygon count and the
complexity of lighting effects etc., while the client only renders the lower
quality version. The pixel-wise colour difference between the versions is then
sent to the client for recovering the higher quality result. The rendering
workload can also be assigned in terms of the number of edges in scene models
(Chen and El-Zarki, 2019) and the number of frames to render (Cuervo et al.,
2015).

Similarly, we adopt a distributed rendering approach but for a hybrid
rendering pipeline. Instead of performing both ray tracing and rasterization
on the server, we offload ray tracing to the server while performing
rasterization locally, combining the results. This also allows for better
utilization of the client hardware as it can still perform rasterization
locally and meet interactive frame rates.

In current cloud gaming services, the gameplay of the user is susceptible to
the stability of the network. Lag and jitter can result which are detrimental
to the interactiveness of metaverse experiences. This is especially crucial
for standalone XR devices used on the go where network connections can be
unstable as compared to in a controlled environment such as their home with a
hardwired Ethernet connection or a fast wireless router. Nonetheless, we
maintain a baseline of temporal coherence for the user by falling back on a
lower quality result if the network connection is unstable. We obtain this
result by predicting the data from the server if it does not arrive on time
and combining it with local rasterization, generating approximate yet up-to-
date interactive output for the user even under poor network conditions. As
such, although our approach works best under fast and stable connections, it
can also handle slow and unstable connections well unlike fully remote
rendering. Cuervo et al. (2015) allows both online and offline gaming.
However, we handle variable network delay and jitter instead of just mid-game
disconnections and provide an approximate prediction of the actual output
instead of simply displaying the local render as-is.

Our technique, Cuervo et al. (2015)’s delta encoding approach, Chen and El-
Zarki (2019) as well as Chuah and Cheung (2014) and Chuah et al. (2016) can
always locally render frames with geometrically accurate scene representation.
The client already has all the data it needs to interactively display a
reasonable low-quality output to the user in real-time. In contrast, for Chan
et al. (2017) and Cuervo et al. (2015)’s client-side I-frame rendering
technique where a subset of the scene geometry and frames respectively are
only rendered at the server, scene representation as shown in the output will
be outdated in the event of network delay and jitter. This same issue also
exists for fully remote rendering where the client is only set up to display
video output rather than perform any sort of rendering locally. We also save
on bandwidth as compared to remote rendering and the distributed rendering
techniques (Cuervo et al., 2015; Chen and El-Zarki, 2019; Chuah and Cheung,
2014; Chuah et al., 2016) as we send variable-size visibility bitmaps based on
the number of lights in the scene from the server as compared to streaming
fixed-size full-colour buffers to the client. We test our approach on ray-
traced shadows, but there are also global illumination approaches in the same
vein (Crassin et al., 2015; Stengel et al., 2021; Magro et al., 2020) where
indirect lighting is computed on servers and direct lighting is calculated on
clients.

##  3\. Design

We perform ray tracing on remote servers, leveraging hardware-accelerated
GPUs. Nonetheless, thin clients like standalone XR devices can still attain
desirable performance with rasterization. Hence, we also leverage their
limited graphics capability for rasterization. DHR can give us not only
interactive frame rates through raster-based rendering on the client, but also
high-quality graphics through ray tracing on the server. We study our
distributed rendering approach on a simple hybrid rendering pipeline (Wyman,
2018) that augments diffuse Lambertian shading (Koppal, 2014) with ray-traced
shadows.

###  3.1. Shading Model

For deferred shading, rasterization is first performed to generate a G-Buffer,
which is a collection of textures containing the data of the nearest fragment
per pixel required for lighting computation. Next, light visibility
information is obtained via ray tracing for every light per pixel. The
G-Buffer and light visibility information are then used to obtain the final
pixel colour $I$ as shown.

(1)
$I=\frac{I_{d}}{\pi}\sum_{i}k_{i}\cdot\textnormal{{saturate(}}\textnormal{\emph{$N\cdot
L_{i}$}}\textnormal{{)}}I_{i}$

$I_{d}/\pi$ refers to the pixel’s diffuse bidirectional reflectance
distribution function with $I_{d}$ as its material diffuse colour. $N$ denotes
its surface normal, and $k_{i}$, $L_{i}$ and $I_{i}$ refer to the relative
visibility, direction and intensity of light $i$ respectively in relation to
the pixel. For every pixel, we trace a shadow ray from its world position to
every light in the scene. If the ray reaches the light without hitting any
other object, the light is deemed visible and $k_{i}=1$, else $k_{i}=0$.

The pixel world positions can be obtained from rasterization or ray casting.
For rasterization, the server can wait for the client to compute the world
positions and send them over or perform rasterization itself to obtain the
world positions earlier. To minimize data transfer, both the client and server
perform rasterization in our current implementation as it is still relatively
fast on both ends. However, while rasterization might have better performance
for simpler scenes, modern ray tracing-accelerated hardware can query primary
visibility with ray casting even faster than through rasterization. This also
streamlines our workflow and avoids repeated computation by omitting
rasterization on the server altogether.

In the original non-distributed ray-traced shadow algorithm, the colour
contributions of lights are computed and accumulated while looping through
every light and tracing rays to it. DHR separates this colour computation and
the ray tracing loop so they can be performed on different hardware for better
parallelism and performance. As such, the ray tracing process only stores the
visibility boolean $k_{i}$ of every light in a visibility bitmap buffer such
that every pixel is assigned a compact bitmask representing all the lights in
the scene. For instance, if light 2 is visible from a particular pixel,
querying the bitmap at its coordinates will produce a bitmask with bit 1 at
position 2. The number of bitmaps and size of per-pixel bitmasks can be tuned
based on the number of lights in the scene.

We show the hybrid rendering output of the shading model in Figure 1 on the
scene The Modern Living Room (CC BY) by Wig42 (Wig42, 2014) (commonly referred
to as Pink Room). The images on the left represent the visualization of the
ray-traced visibility buffer encapsulating information for 3 scene lights
where the material colour of each pixel is shown if the respective light is
not obstructed from its world space position. The right image is the final
result of combining the ray-traced shadow information in the visibility buffer
with diffuse Lambertian shading at the client.

Figure 1. Shading output.

3 images on the left representing each of the 3 lights display visibility
information while 1 image on the right shows the final output.

###  3.2. Performance-Accuracy Tradeoff

Our DHR approach in Figure 2 adopts the User Datagram Protocol (UDP) for fast
data transmission. Although UDP does not retransmit dropped packets,
retransmission delays tend to be too costly for interactive applications
anyway. Before any rendering is done at the server, the user’s scene inputs
need to be sent from the client to the server for synchronization. For now, we
only handle user-controlled camera movement so the scene inputs consist of
camera data (i.e. position, target and up vector). After computing the
visibility buffer corresponding to this camera data, the server then
compresses it with the LZ4 algorithm which is lossless and efficient on
bitmaps. It then sends the compressed buffer and its corresponding frame
number to the client for shading computation.

Figure 2. Control flow diagram.

The DHR pipeline includes the G-Buffer, prediction, visibility, shading,
network and memory transfer passes.

Taking into account network latency, while the client updates the server with
the camera data for every frame, the frame it displays to the user is
dependent on the latest visibility buffer received from the server. We define
response time as the elapsed time from the client sending the camera data to
receiving its corresponding visibility buffer. As such, although the client
can compute the most updated G-Buffer locally, its latest received visibility
buffer will be around $p$ frames behind, where $p$ is the response time
divided by the total time taken to render and display one frame. For example,
at 90 FPS and a response time of 55 ms (50 ms of network ping and 5 ms of
rendering on the server), the client will have access to a visibility buffer
of approximately $55/11.1{\approx}5$ frames behind its G-Buffer. The wait for
the correct updated visibility buffer might take too long if the network is
slow or even be indefinite if the data is lost as there will be no
retransmission from the server.

To strike a balance between performance and accuracy, we provide some
allowance for the visibility buffer used such that it can be behind the
G-Buffer by a user-defined maximum number of frames $x_{max}$. The value of
$x_{max}$ is an artistic choice set by the game developer based on the speed
of objects and camera movement in the scene. This is a reasonable
approximation as even for dynamic scenes, frames temporally near to each other
should be spatially similar, in this case in terms of the occlusion of
geometry with respect to scene lights. Nonetheless, applying an outdated
visibility buffer as-is with a more recent G-Buffer would lead to shadows not
perfectly matching the scene geometry in the final image.

###  3.3. Prediction

Hence, to improve the accuracy of the visibility buffer used, we perform a
prediction of the correct visibility buffer based on the scene information
with respect to the G-Buffer and offset the visibility buffer to match this
scene. In doing so, we minimize the error resulting from scene misalignment
between the two buffers.

In Figure 2, the client is on its $n^{\text{th}}$ frame with a latest
visibility buffer of the $(n-p)^{\text{th}}$ frame, so it can predict
$x=min(x_{max},p)$ frames to obtain an approximation of the
$(n-p+x)^{\text{th}}$ visibility buffer and render the more recent
$(n-p+x)^{\text{th}}$ frame with the $(n-p+x)^{\text{th}}$ G-Buffer. This
allows the client to experience a faster ”response time” while keeping the
amount of approximation to a pre-defined minimum, maintaining a consistent
level of visual quality for the user. For instance, setting $x=2$ in the
example from Subsection 3.2 gives a small lag of around 33 ms or 5 - 2 = 3
frames.

For the prediction, the client writes the camera’s view projection matrix for
every frame in a circular buffer which can store up to the buffer size amount
of the latest consecutive frames. As shown in Figure 3, it takes the latest
visibility buffer along with its corresponding frame number to perform frame
prediction. The prediction pass uses the received frame number to retrieve the
older camera view projection matrix corresponding to the received visibility
buffer. Together with the newer camera view projection matrix of the G-Buffer,
the prediction pass calculates a two-dimensional motion vector for every
pixel. These motion vectors store the difference between the screen space
position of the same world space points as seen by the old camera and the new
camera. Finally, each pixel’s visibility information in the received
visibility buffer will be offset by the individual motion vectors, forming a
new visibility buffer that will be used for the shading of the scene.

Figure 3. Visibility buffer prediction pipeline.

If the prediction is executed, world space data is used to generate screen-
space motion vectors. The motion vectors are then used to offset the received
visibility buffer.

We do not know the visibility information of scene points revealed by camera
movement, so their corresponding pixels are taken to be fully illuminated in
relation to all scene lights. Doing so helps to prevent adding false shadows
to the scene if there is camera movement. However, there are prominent borders
of illuminated pixels along the edges of the output frame in the direction of
camera movement as the offset positions given by the motion vectors fall
outside the boundary of the original visibility buffer. We minimize this error
at a small performance cost by working with slightly larger textures so some
visibility information is available past the borders of the frame. For
example, taking a display resolution of 1920 $\times$ 1080, we work with a
larger 1984 $\times$ 1116 visibility buffer texture with an additional 64
$\times$ 36 pixels. This visibility buffer size is another variable that the
game developer can set depending on the interactive performance requirement of
the application and the speed of camera movement. The focal length of the
camera is adjusted to maintain the same displayed region of the scene for
different numbers of additional pixels to the display resolution.

##  4\. Evaluation

The results of our DHR implementation are measured on a client as well as a
server with the ray tracing-accelerated NVIDIA GeForce RTX 2080 GPU. Both the
client and server have an Intel Core i7-7700K central processing unit (CPU)
and 16 GB random-access memory (RAM). Our qualitative measurements are taken
over a 14 second and 50 second camera animation sequence for the scenes Pink
Room (PR) and UE4 Sun Temple (CC BY-NC-SA) by Epic Games (2017) (ST)
respectively at a display resolution of 1920 $\times$ 1080.

###  4.1. Graphics Quality

We show the result of our visibility buffer prediction on the final output
with camera movement. As seen in Table 1, our prediction pass improves the
accuracy of the shadows. With an additional 128 $\times$ 72 pixels on the
visibility buffer, we obtain Figures 4(d) and 5(d) and their corresponding
objective visual metrics measurements as shown in Figure 6. The visual quality
of the output decreases as more frames $x$ are predicted for the same frame
$n$. This is because as $x$ increases while $n$ remains constant, $n-p$
decreases which means that an older visibility buffer is used for the
prediction, resulting in a higher bitwise prediction error as shown in Figure
7. From our preliminary testing with some users, $x$ is best capped at 4 or 5.

Table 1. Prediction output. Delay (ms) | 0 | 100 | 200
---|---|---|---
Before Prediction |  |  |
After Prediction |  |  |

(a) No prediction (b) $x=1$

(c) $x=5$ (d) $x=10$ Figure 4. PR (454,505 vertices).

Pink Room, medium shot of coffee table with sofa in background

(a) No prediction (b) $x=1$

(c) $x=5$ (d) $x=10$ Figure 5. ST (1,641,711 vertices).

Sun Temple, medium shot of angel casting a shadow on the background wall

SSIM, PSNR and VMAF values decrease as the number of frames predicted
increases.

Figure 6. Objective visual metrics.

Figure 7. Prediction error.

Error is large when the visibility buffer is small and the number of frames
predicted is large.

Under the fixed laboratory testing environment used for our experiments, the
Structural Similarity Index Measure (SSIM), Peak Signal-to-Noise Ratio (PSNR)
and Video Multimethod Assessment Fusion (VMAF) scores obtained do not vary
across repeated experiments. The VMAF data also does not highlight the impact
of the visibility buffer size on the accuracy of the prediction. As such,
there is room for improvement in finding a better metric to measure the human
perception of error arising from the misalignment of the visibility buffer.
Nonetheless, there is a general downward trend at a low number of frames
predicted, showing that human perception of the prediction errors can be felt.
However, there are anomalies when more frames are predicted. This could be
because VMAF is mostly trained to identify and evaluate the effects of
artifacts and blurriness arising from lossy encoding. While the artifacts
caused by the prediction error may resemble such artifacts when small, VMAF
may not be able to accurately evaluate the video when the prediction error and
hence image distortion becomes significant. Regardless, even amongst the
anomalies, we can observe a downward trend.

Figure 7 shows the average prediction error per light measured using the
bitwise comparison of the actual visibility buffer with the predicted version.
We see that increasing the size of the visibility buffer does lower the
prediction error by roughly a constant amount, which can only arise from
minimizing errors at the edges of the frame since increasing buffer size does
not affect anything else. This has a large impact and is visually very
noticeable at lower values of frames predicted, where the error can be almost
halved with the use of 128 $\times$ 72 additional pixels. We included a result
of using an unreasonably large buffer of 2560 $\times$ 1440 pixels (i.e. 640
$\times$ 360 additional pixels) on PR to show that a large portion of the
error arises from the edges of the frame. The same result was not measured on
the larger ST due to the difficulty of maintaining the frame rate while
performing the error calculation.

The prediction error calculations used samples recorded at a fixed frame rate
of 60 FPS as an industry standard frame rate for real-time interactive
applications like games. VMAF samples were recorded at a lower 30 FPS to
reduce deviation from VMAF’s 24 FPS frame rate recommendation. A discrepancy
from the recommended value would not affect the reading significantly, but it
will be interpreted as a 24 FPS video so the motion feature of the metric will
generate lower scores and cause the overall VMAF score to be slightly lower.

Qualitatively, the shadows for both scene setups generally remain consistent
except for those cast by foreground geometry onto background objects. For
example, the moving shadows of the red vase in Figure 8(b) and the angel
statue in Figure 9(b) with different numbers of frames predicted show how the
visibility information of fragments in the interior of the frame revealed by
camera movement cannot be reliably predicted by offsetting the visibility
buffer of previous frames which do not contain this information in the first
place. As for pixels at the edges of the frame, Figure 4(d) does not have the
issue of illuminated pixel borders like Figure 5(d), implying that the ST
scene setup here needs to be rendered with a larger visibility buffer texture.
Additionally, in Figure 8(b), we see jaggies along the shadow edge as it is
enlarged from the original visibility buffer to the offset one due to camera
movement, so it now occupies more pixels. This results in multiple nearby
pixels querying the original visibility buffer at the same coordinates which
can be avoided with linear sampling instead of point sampling. We can also see
that in Figure 4(d), the shadow on the floor of the beige-coloured mug which
is moving in PR is not predicted accurately as we currently do not handle
dynamic scene geometry in our visibility buffer prediction.

(a) $x=1$ (b) $x=10$ Figure 8. Figure 4(d) (PR) zoomed in.

Pink Room, close-up shot of red vase on coffee table

(a) $x=1$ (b) $x=10$ Figure 9. Figure 5(d) (ST) zoomed in.

Sun Temple, close-up shot of right angel wing

###  4.2. Performance

In our setup, the client and server are connected via Ethernet over a local
network and are physically located very close to each other (i.e. in the same
lab). As such, the latency between them is less than 1 ms. Hence, we add
artificial network ping via network condition simulator _clumsy_ and measure
the frame rate and response time of our approach on different network latency
conditions. The results in Table 2 are as expected as they show that frame
rate decreases with increased geometric complexity of the scene and additional
pixels of the visibility buffer. In general, we achieve interactive frame
rates for our scene setups even when rendering large visibility buffers.

Table 2. Frame rate. Additional Pixels | PR | ST
---|---|---
0 | 154 | 133
64 $\times$ 36 | 148 | 127
128 $\times$ 72 | 142 | 125

As shown in Figure 10, the response time is fast for low network ping values.
Initially at no frames predicted, the response time is higher with a larger
buffer as it takes slightly longer for the server to transfer between device
and host memory, render, compress and send the texture. Response time then
decreases as expected with more frames predicted. We also note that the
response time experienced when using different visibility buffer sizes
converges as seen in the bottom graph. This is because every additional frame
predicted reduces the response time more for larger visibility buffers which
take more time to generate and send to the client.

Figure 10. Response time.

Response time decreases as the number of frames predicted increases

##  5\. Future Work

To improve the accuracy of our output, one possible modification to the
pipeline could be to perform server-side prediction (Anand and Wenren, 2017)
in addition to the current client-side prediction. The server would render a
few different visibility buffers based on cases of potential camera movement.
These visibility buffers would be almost identical and contain much
duplication of data, so they should be able to be highly compressed before
they are sent to the client. The client will then select the visibility buffer
that corresponds most closely to the actual camera motion and use that for the
client-side prediction.

Currently, if the server cannot be reached at all, the client erroneously
performs lighting computation with an empty visibility buffer. However, in
extreme cases of poor network conditions, we are intending to fall back on
local rendering with just rasterization as a last resort. This enables us to
generate interactive output for the user even if the connection is too slow or
lost, or if the remote server is down. We would be able to maintain a baseline
of temporal coherence for the user in such situations by switching to local
rendering dynamically when network latency values exceed a certain threshold,
which could be a variable we expose to developers for trading-off between
visual quality and performance.

Performing rasterization on the client while ray tracing on the server in
parallel can give an overall improvement in performance as compared to pure
remote rendering when enabled by fast 5G networks and low latencies to reach
5G edge nodes (Anand and Hao Edwin, 2014). Hence, we intend to test our
implementation in a simulated 5G edge computing environment to leverage the
fast network technology.

We are also looking to expand the scope of scene inputs to accommodate dynamic
changes in scene geometry in addition to camera movement. Besides affecting
the scene synchronization between the client and server, this upgrade will
also impact the prediction pass. On top of the camera’s view projection
matrix, the circular buffer would also need to hold the information of moving
objects for every frame, which would have to be taken into account when
calculating the screen space motion vector of each pixel.

The compression quality of visibility buffer packets of the same frame can
also be varied to save on bandwidth. The user’s attention is more likely to
focus on certain regions in the frame so Babaei et al. (2017) maps out
predictions of such regions and allocates more bitrate accordingly, while
Illahi et al. (2020) employs a gaze tracker device to determine the exact spot
the user is looking at.

In the future, we hope to incorporate multiple complex lighting and camera
effects into a single hybrid rendering pipeline as part of an adaptive hybrid
system that can dynamically adjust parameters for the different effects to
achieve the maximum visual quality while maintaining interactive frame rates,
being adaptive to changing scene properties and hardware configurations. We
also hope for this pipeline to adopt a DHR approach for ray tracing-
incorporated graphics in standalone XR devices as part of metaverse
experiences.

##  6\. Conclusion

We present the Distributed Hybrid Rendering (DHR) approach for standalone
extended reality (XR) devices through the design and evaluation of a simple
prototype which generates ray-traced shadows without compromising the
interactive frame rates required for immersive metaverse experiences. Our
technique is elastic to network conditions, dynamically adjusting the amount
of visibility buffer approximation to maximize the visual quality of the
output while maintaining interactive frame rates.

###### Acknowledgements.

We thank Low Siang Ern for helping to improve the writing and quality of the
paper. This work is supported by the Singapore Ministry of Education Academic
Research grant T1 251RES2205, “Real-time Distributed Hybrid Rendering with 5G
Edge Computing for Realistic Graphics in Mobile Games and Metaverse
Applications”.

## References

  * (1)
  * Amazon (2021) Amazon. 2021.  Luna.  https://www.amazon.com/luna/
  * Anand and Hao Edwin (2014) Bhojan Anand and Aw Jia Hao Edwin. 2014.  Gamelets — Multiplayer mobile games with distributed micro-clouds. In _2014 Seventh International Conference on Mobile Computing and Ubiquitous Networking (ICMU)_. 14–20.  https://doi.org/10.1109/ICMU.2014.6799051
  * Anand and Wenren (2017) Bhojan Anand and Pan Wenren. 2017.  CloudHide: Towards Latency Hiding Techniques for Thin-Client Cloud Gaming. In _Proceedings of the on Thematic Workshops of ACM Multimedia 2017_ (Mountain View, California, USA) _(Thematic Workshops ’17)_. Association for Computing Machinery, New York, NY, USA, 144–152.  https://doi.org/10.1145/3126686.3126777
  * Babaei et al. (2017) Ebrahim Babaei, Mahmoud Reza Hashemi, and Shervin Shirmohammadi. 2017.  A State-Based Game Attention Model for Cloud Gaming. In _Proceedings of the 15th Annual Workshop on Network and Systems Support for Games_ (Taipei, Taiwan) _(NetGames ’17)_. IEEE Press, 34–36.
  * Barré-Brisebois et al. (2019) Colin Barré-Brisebois, Henrik Halén, Graham Wihlidal, Andrew Lauritzen, Jasper Bekkers, Tomasz Stachowiak, and Johan Andersson. 2019.  Hybrid Rendering for Real-Time Ray Tracing.  In _Ray Tracing Gems_ , Eric Haines and Tomas Akenine-Möller (Eds.). Apress, Chapter 25.  http://raytracinggems.com.
  * Beck et al. (1981) Steve Beck, Andreas-Christoph Bernstein, Daniel Danch, and Bernd Fröhlich. 1981\.  CPU-GPU Hybrid Real Time Ray Tracing Framework, Vol. 0. The Eurographics Association and Blackwell Publishing Ltd., 1–8.  https://api.semanticscholar.org/CorpusID:18758412
  * Cabeleira (2010) João Pedro Guerreiro Cabeleira. 2010\.  _Combining Rasterization and Ray Tracing Techniques to Approximate Global Illumination in Real-Time_.  Master’s thesis. Portugal.  http://voltaico.net/files/article.pdf
  * Chan et al. (2017) Kar-Long Chan, K. Ichikawa, Yasuhiro Watashiba, Putchong Uthayopas, and Hajimu Iida. 2017\.  A Hybrid-Streaming Method for Cloud Gaming: To Improve the Graphics Quality delivered on Highly Accessible Game Contents.  _Int. J. Serious Games_ 4 (2017).  https://doi.org/10.17083/ijsg.v4i2.163
  * Chen and El-Zarki (2019) De-Yu Chen and Magda El-Zarki. 2019.  A Framework for Adaptive Residual Streaming for Single-Player Cloud Gaming.  _ACM Trans. Multimedia Comput. Commun. Appl._ 15, 2s, Article 66 (July 2019), 23 pages.  https://doi.org/10.1145/3336498
  * Chuah and Cheung (2014) Seong-Ping Chuah and Ngai-Man Cheung. 2014.  Layered Coding for Mobile Cloud Gaming. In _Proceedings of International Workshop on Massively Multiuser Virtual Environments_. Association for Computing Machinery, New York, NY, USA, 1–6.  https://dl.acm.org/doi/10.1145/2577387.2577395
  * Chuah et al. (2016) Seong-Ping Chuah, Ngai-Man Cheung, and Chau Yuen. 2016.  Low Bit-Rate Mobile Cloud Gaming. In _Proceedings of the 3rd Workshop on Mobile Gaming_ (Singapore, Singapore) _(MobiGames ’16)_. Association for Computing Machinery, New York, NY, USA, 17–22.  https://doi.org/10.1145/2934646.2934649
  * Crassin et al. (2015) Cyril Crassin, David Luebke, Michael Mara, Morgan McGuire, Brent Oster, Peter Shirley, Peter-Pike Sloan, and Chris Wyman. 2015\.  CloudLight: A System for Amortizing Indirect Lighting in Real-Time Rendering.  _Journal of Computer Graphics Techniques (JCGT)_ 4, 4 (15 October 2015), 1–27.  http://jcgt.org/published/0004/04/01/
  * Cuervo et al. (2015) Eduardo Cuervo, Alec Wolman, Landon Cox, Kiron Lebeck, Ali Razeen, Stefan Saroiu, and Madanlal Musuvathi. 2015.  Kahawai: High-Quality Mobile Gaming Using GPU Offload. In _Proceedings of the 13th Annual International Conference on Mobile Systems, Applications, and Services_ (Florence, Italy) _(MobiSys ’15)_. Association for Computing Machinery, New York, NY, USA, 121–135.  https://doi.org/10.1145/2742647.2742657
  * Eisert and Fechteler (2008) Peter Eisert and Philipp Fechteler. 2008.  Low Delay Streaming of Computer Graphics. In _2008 15th IEEE International Conference on Image Processing_. 2704–2707.  https://doi.org/10.1109/ICIP.2008.4712352
  * Epic Games (2017) Epic Games. 2017\.  Unreal Engine Sun Temple, Open Research Content Archive (ORCA).  https://developer.nvidia.com/ue4-sun-temple
  * Google (2021) Google. 2021.  Stadia.  https://stadia.google.com/
  * Hertel et al. (2009) Stefan Hertel, Kai Hormann, and Rüdiger Westermann. 2009\.  A Hybrid GPU Rendering Pipeline for Alias-Free Hard Shadows. In _Eurographics 2009 Areas Papers_ , D. Ebert and J. Krüger (Eds.). München, Germany, 59–66.
  * Holthe et al. (2009) Ole-Ivar Holthe, Ola Mogstad, and Leif Arne Ronningen. 2009\.  Geelix LiveGames: Remote Playing of Video Games. In _2009 6th IEEE Consumer Communications and Networking Conference_. 1–2.  https://doi.org/10.1109/CCNC.2009.4784713
  * Illahi et al. (2020) Gazi Karam Illahi, Thomas Van Gemert, Matti Siekkinen, Enrico Masala, Antti Oulasvirta, and Antti Ylä-Jääski. 2020\.  Cloud Gaming with Foveated Video Encoding.  _ACM Trans. Multimedia Comput. Commun. Appl._ 16, 1, Article 7 (Feb. 2020), 24 pages.  https://doi.org/10.1145/3369110
  * Koppal (2014) Sanjeev Koppal. 2014\.  _Lambertian Reflectance_.  Springer US, Boston, MA, 441–443.  https://doi.org/10.1007/978-0-387-31439-6_534
  * Lauterbach and Manocha (2009) Christian Lauterbach and Dinesh Manocha. 2009.  _Fast Hard and Soft Shadow Generation on Complex Models Using Selective Ray Tracing_.  Technical Report TR09-004. UNC CS.
  * Lee et al. (2021) Lik-Hang Lee, Tristan Braud, Pengyuan Zhou, Lin Wang, Dianlei Xu, Zijun Lin, Abhishek Kumar, Carlos Bermejo, and Pan Hui. 2021.  All One Needs to Know about Metaverse: A Complete Survey on Technological Singularity, Virtual Ecosystem, and Research Agenda.  https://doi.org/10.48550/ARXIV.2110.05352
  * Macedo et al. (2018) Daniel Valente De Macedo, Ygor Rebouças Serpa, and Maria Andréia Formico Rodrigues. 2018.  Fast and Realistic Reflections Using Screen Space and GPU Ray Tracing—A Case Study on Rigid and Deformable Body Simulations.  _Comput. Entertain._ 16, 4, Article 5 (Nov. 2018), 18 pages.  https://doi.org/10.1145/3276324
  * Magro et al. (2020) Mark Magro, Keith Bugeja, Sandro Spina, and Kurt Debattista. 2020\.  Cloud-Based Dynamic GI for Shared VR Experiences.  _IEEE Computer Graphics and Applications_ 40, 5 (2020), 10–25.  https://doi.org/10.1109/MCG.2020.3002474
  * Marrs et al. (2018) Adam Marrs, Josef Spjut, Holger Gruen, Rahul Sathe, and Morgan McGuire. 2018.  Adaptive Temporal Antialiasing. In _Proceedings of the Conference on High-Performance Graphics_ (Vancouver, British Columbia, Canada) _(HPG ’18)_. ACM, New York, NY, USA, Article 1, 4 pages.  https://doi.org/10.1145/3231578.3231579
  * NVIDIA (2021) NVIDIA. 2021.  GeForce NOW.  https://www.nvidia.com/en-us/geforce-now/
  * Parisi (2021) Tony Parisi. 2021\.  The Seven Rules of the Metaverse.  https://medium.com/meta-verses/the-seven-rules-of-the-metaverse-7d4e06fa864c
  * Stengel et al. (2021) Michael Stengel, Zander Majercik, Benjamin Boudaoud, and Morgan McGuire. 2021.  A Distributed, Decoupled System for Losslessly Streaming Dynamic Light Probes to Thin Clients. In _Proceedings of the 12th ACM Multimedia Systems Conference_ (Istanbul, Turkey) _(MMSys ’21)_. Association for Computing Machinery, New York, NY, USA, 159–172.  https://doi.org/10.1145/3458305.3463379
  * Tan et al. (2021) Yu Wei Tan, Louiz Kim-Chan, Anthony Halim, and Anand Bhojan. 2021.  Cloud-Assisted Hybrid Rendering for Thin-Client Games and VR Applications. In _Pacific Graphics Short Papers, Posters, and Work-in-Progress Papers_ , Sung-Hee Lee, Stefanie Zollmann, Makoto Okabe, and Burkhard Wünsche (Eds.). The Eurographics Association.  https://doi.org/10.2312/pg.20211389
  * Wig42 (2014) Wig42. 2014.  The Modern Living Room.  https://www.blendswap.com/blends/view/75692
  * Wyman (2018) Chris Wyman. 2018\.  Introduction to DirectX Raytracing. In _ACM SIGGRAPH 2018 Courses_ (Vancouver, Canada) _(SIGGRAPH ’18)_.  http://cwyman.org/code/dxrTutors/dxr_tutors.md.html


11institutetext:  University of Manchester, Manchester, UK 22institutetext:
Institute of Public Administration, Jeddah, Saudi Arabia

# FuSeBMC v4: Smart Seed Generation for Hybrid Fuzzing

Kaled M. Alshmrany 1122 Mohannad Aldughaim 11 Ahmed Bhayat 11
Lucas C. Cordeiro 11

###### Abstract

FuSeBMC is a test generator for finding security vulnerabilities in C
programs. In earlier work [4], we described a previous version that
incrementally injected labels to guide Bounded Model Checking (BMC) and
Evolutionary Fuzzing engines to produce test cases for code coverage and bug
finding. This paper introduces a new version of FuSeBMC that utilizes both
engines to produce smart seeds. First, the engines are run with a short time
limit on a lightly instrumented version of the program to produce the seeds.
The BMC engine is particularly useful in producing seeds that can pass through
complex mathematical guards. Then, FuSeBMC runs its engines with more extended
time limits using the smart seeds created in the previous round. FuSeBMC
manages this process in two main ways using its Tracer subsystem. Firstly, it
uses _shared memory_ to record the labels covered by each test case. Secondly,
it evaluates test cases, and those of high impact are turned into seeds for
subsequent test fuzzing. As a result, we significantly increased our code
coverage score from last year, outperforming all tools that participated in
this year’s competition in every single category.

##  1 Overview

Software testing is one of the most crucial phases in software development
[11]. Tests often expose critical bugs in software applications. In earlier
work [4], we presented FuSeBMC, an automated test generation tool that
exploits the combination of Fuzzing and Bounded Model Checking. FuSeBMC
achieved second place in Test-Comp 2021 [5, 3] and first place in the Cover-
Error category. It ranked fourth in the Cover-Branches category. This year, we
introduce a new version of FuSeBMC (v4) that adds smart seed generation and
shared memory amongst other improvements and features. The new version
significantly improves on the previous version, particularly relating to code
coverage. One of the primary contributions of this paper is the linking of a
grey-box fuzzer with a bounded model checker. A bounded model checker works by
treating a program as a state transition system and then checking whether
there exists a transition in this system of length less than a bound $k$ that
violates the property to be verified [6, 8]. We leverage this power of model
checkers as a method for smart seed generation. During grey-box fuzzing, if a
particular branch has not been explored, BMC can be used to provide a model
(set of assignments to input variables), which reaches the branch. This model
is then shared in memory and seeded for further grey-box fuzzing. However, BMC
can be slow and resource-intensive. We also carry out a lightweight static
program analysis to recognize input verification. We analyze the code for
conditions on the input variables and ensure that seeds are only selected if
they pass these conditions. Together, these contributions turn FuSeBMC into a
world-leading fuzzer.

##  2 Test Generation Approach

Figure 1 provides an overview of the components within FuSeBMC and how these
interact. FuSeBMC makes use of the Clang tooling infrastructure [1] to
instrument programs. In addition, FuSeBMC employs three engines in its
reachability analysis. FuSeBMC starts by injecting goal labels into the given
input (C program) and ranks them depending on the chosen strategy. After that,
FuSeBMC employs two fuzzers, one based on the American Fuzzy Loop (AFL) [7,
2], and a second custom fuzzer, which we refer to as _selective fuzzer_ (see
[4] for details). ESBMC [10, 9] is a state-of-the-art SMT-based bounded model
checker which FuSeBMC utilises to produce seeds and test cases. Besides for
these engines, FuSeBMC also incorporates a subsystem we call the Tracer. This
component maintains a program graph and records which labels in the graph have
been covered by which test cases. It coordinates the other tools and handles
the passing of information between them. The fuzzers generate test cases by
randomly mutating the program’s input and running it to analyze code coverage.
The BMC engine executes the given program symbolically to determine the
reachability of particular goal labels. In case of success, a witness – a set
of inputs leading to the goal – is produced. If any engine manages to reach a
label, we say that label is _covered_. The Tracer records which goals have
been covered by which test cases. This information is used to prevent the
computationally expensive BMC engine from trying to reach an already covered
goal. More importantly, it is used to coordinate the fuzzing and BMC engines.
Let $L_{n}$ be the deepest label on some branch that ESBMC covers and let
$L_{n+1}$ be the next deepest label on the branch. Tracer records the test
case produced by ESBMC that covers $L_{n}$ and then passes this as a seed to
the fuzzer. Mutating such a seed has a far larger chance of leading to a test
case that covers $L_{n+1}$ than starting the fuzzer from scratch. FuSeBMC runs
until all goals are covered or a timeout is reached.

Figure 1: FuSeBMC v4 Framework. This figure displays the major components of
the FuSeBMC test generator and how they interact. Note in particular the seed
store.

###  2.0.1 Code Instrumentation

FuSeBMC front-end uses Clang tooling infrastructure [1] to parse a C program
and produce an Abstract Syntax Tree (AST). While traversing the AST, FuSeBMC
injects labels into each branch, including every conditional statement, loop,
and function. Using these labels, FuSeBMC can measure the code coverage.

###  2.0.2 Reachability Graph Analysis

After instrumenting the C program, FuSeBMC analyzes it and produces a
reachability graph. The graph will assign each goal label to the code block it
is located in. Then, FuSeBMC ranks goals depending on the strategy chosen. For
example, one strategy is to prefer deeper goals over shallower goals. This
strategy improves the performance of FuSeBMC since a test case that covers a
deep goal will also cover shallower goals on the path to it. FuSeBMC also
ranks coverage metrics over others, such as conditional coverage over loop
coverage.

###  2.0.3 Seed Generation

A unique aspect of the latest version of FuSeBMC is a seed generation phase
that is run prior to the start of the principal reachability analysis. In this
phase, FuSeBMC first lightly instruments the code under test by limiting loop
bounds and assuming a narrow range of values for input variables. The bounds
on input variables are further limited by carrying out a lightweight static
analysis to recognize code that applies verification conditions to input
variables. After instrumenting the code, FuSeBMC runs its fuzzing and BMC
engines with very short time limits. The test cases generated by these engines
are ranked, and the highest impact test cases are selected as smart seeds for
the next round. The impact of a test case is measured using two metrics.
First, the number of labels covered uniquely by that test case, and second,
the program depth achieved by the test case. ESBMC is particularly effective
at seed generation as its underlying SMT solvers can be used to discover test
cases that circumvent complex mathematical guards.

###  2.0.4 Reachability Analysis Engines

In its primary phase, FuSeBMC carries out reachability analysis. Essentially
this involves running the engines with longer timeouts on the original non-
instrumented code with the fuzzer making use of the smart seeds. The Tracer
coordinates the various engines through the use of _shared memory_. For
example, assume that ESBMC is unable to cover some goal $L$ and let
$L^{\prime}$ be the deepest goal on the path to $L$ that ESBMC _can_ cover.
The tracer records the test case that covers $L^{\prime}$ and passes it to the
fuzzer via shared memory as an incomplete seed. Thus, FuSeBMC combines the
strengths of both types of engines. The BMC engine produces seeds that bypass
complex guards and thereby help the fuzzers explore paths deep within the
program. During reachability analysis, the Tracer constantly monitors the test
cases produced by the various engines. High impact test cases, as measured by
the metrics discussed above, are passed by the Tracer back to the seed store.
Thus, the seed store is dynamically updated as the analysis progresses.

##  3 Strengths and Weaknesses

The strengths of the latest version of FuSeBMC are as follows. It runs a
dedicated seed generation phase in order to start the main fuzzing effort with
high quality, high impact seeds. Furthermore, during the main test-generation
phase, these seeds are constantly being updated. Beyond this, it incorporates
a dedicated subsystem, the Tracer, that uses a shared memory store to manage
the various engines. By combining the engines, the Tracer ensures that FuSeBMC
far outperforms the individual engines, or even the running of the engines in
parallel, but isolated. The outcome of these improvements can be seen in the
ECA and Combination benchmark sets. Previously, these posed a challenge to
FuSeBMC. With the latest changes, according to preliminary results, FuSeBMC
achieved first place in the Combination subcategory and took second place in
the ECA subcategory of the 2022 Test-Comp competition. In the ECA subcategory
FuSeBMC displayed a remarkable $60$% improvement over its performance from the
previous competition. The preliminary results also indicate that FuSeBMC has
achieved first place in the Cover-Branches category with high coverage and
validation statistics. One of the weaknesses of FuSeBMC that we plan to work
on, is that for large programs, particularly for programs that redefine C
library functions, seed generation can be slow and consume too much of the
tool’s time.

##  4 Tool Setup and Configuration

FuSeBMC can be run using the command below. The user is required to set the
architecture, the property file path, the competition strategy, and the
benchmark path, as:

fusebmc.py [-a {32, 64}] [-p PROPERTY_FILE]
this text is inv [-s {kinduction,falsi,incr,fixed}]
this text is inv [BENCHMARK_PATH]

where -a sets the architecture to 32 or 64, -p sets the property file to
PROPERTY_FILE, where it has a list of all the properties to be tested. -s sets
the BMC strategy to one of the listed strategies{kinduction,falsi,incr,fixed}.
The Benchexec tool info module is fusebmc.py and the benchmark definition file
is FuSeBMC.xml.

##  5 Software Project

FuSeBMC is implemented using C++, and it is publicly available under the terms
of the MIT License at GitHub111https://github.com/kaled-alshmrany/FuSeBMC. The
repository includes the latest version of FuSeBMC (version 4.1.14). FuSeBMC
dependencies and instructions for building from source code are all listed in
the README.md file.

## References

  * [1] Clang documentation.  http://clang.llvm.org/docs/index.html, 2015.
  * [2] American fuzzy lop, https://lcamtuf.coredump.cx/afl/, 2021.
  * [3] Kaled Alshmrany et al.  FuSeBMC: A white-box fuzzer for finding security vulnerabilities in c programs.  (FASE), 12649:363–367, 2020.
  * [4] Kaled Alshmrany et al.  FuSeBMC: An energy-efficient test generator for finding security vulnerabilities in c programs.  In International Conference on TAP, pages 85–105, 2021.
  * [5] Dirk Beyer.  Status report on software testing: Test-comp 2021.  Proc. FASE. LNCS, 12649.
  * [6] Armin Biere.  Bounded model checking.  Frontiers in Artificial Intelligence and Applications, pages 457–481. 2009.
  * [7] Marcel Böhme et al.  Directed greybox fuzzing.  In Proceedings of the 2017 ACM SIGSAC Conference on Computer and Communications Security, 2017.
  * [8] Lucas C. Cordeiro, Bernd Fischer, and João Marques-Silva.  Smt-based bounded model checking for embedded ANSI-C software.  IEEE Trans. Software Eng., 38(4):957–974, 2012.
  * [9] Mikhail R Gadelha et al.  ESBMC: Scalable and precise test generation based on the floating-point theory:(competition contribution).  FASE 2020, pages 525–529.
  * [10] Mikhail R Gadelha et al.  ESBMC v6. 0: Verifying c programs using k-induction and invariant inference.  In International Conference on TACAS, pages 209–213. Springer, 2019\.
  * [11] Nicha Kosindrdecha and Jirapun Daengdej.  A test case generation process and technique.  journal of Software Engineering, 4(4):265–287, 2010.