# Anisotropic Optical And Vibrational Properties Of GeS

Natalia Zawadzka  Institute of Experimental Physics, Faculty of Physics,
University of Warsaw, ul. Pasteura 5, 02-093 Warsaw, Poland  Łucja Kipczak
Institute of Experimental Physics, Faculty of Physics, University of Warsaw,
ul. Pasteura 5, 02-093 Warsaw, Poland  Tomasz Woźniak  Department of
Semiconductor Materials Engineering, Wrocław University of Science and
Technology, ul. Wybrzeże Wyspiańskiego 27, 50-370 Wrocław, Poland  Katarzyna
Olkowska-Pucko  Institute of Experimental Physics, Faculty of Physics,
University of Warsaw, ul. Pasteura 5, 02-093 Warsaw, Poland  Magdalena
Grzeszczyk  Institute of Experimental Physics, Faculty of Physics, University
of Warsaw, ul. Pasteura 5, 02-093 Warsaw, Poland  Adam Babiński  Institute of
Experimental Physics, Faculty of Physics, University of Warsaw, ul. Pasteura
5, 02-093 Warsaw, Poland  Maciej R. Molas<EMAIL_ADDRESS>Institute
of Experimental Physics, Faculty of Physics, University of Warsaw, ul.
Pasteura 5, 02-093 Warsaw, Poland

###### Abstract

The optical response of bulk germanium sulfide (GeS) is investigated
systematically using different polarization-resolved experimental techniques,
such as photoluminescence (PL), reflectance contrast (RC), and Raman
scattering (RS). It is shown that while the low-temperature ($T$=5 K) optical
band-gap absorption is governed by a single resonance related to the neutral
exciton, the corresponding emission is dominated by the disorder/impurity-
and/or phonon-assisted recombination processes. Both the RC and PL spectra are
found to be linearly polarized along the armchair direction. The low and room
($T$=300 K) temperature RS spectra consist of six Raman peaks identified with
the help of Density Fuctional Theory (DFT) calculations:
A${}^{1}_{\textrm{g}}$, A${}^{2}_{\textrm{g}}$, A${}^{3}_{\textrm{g}}$,
A${}^{4}_{\textrm{g}}$, B${}^{1}_{\textrm{1g}}$, and B${}^{2}_{\textrm{1g}}$,
which polarization properties are studied under four different excitation
energies. We found that the polarization orientations of the
A${}^{2}_{\textrm{g}}$ and A${}^{4}_{\textrm{g}}$ modes under specific
excitation energy can be useful tools to determine the GeS crystallographic
directions: armchair and zigzag.

##  I Introduction

Two-dimensional (2D) layered van der Waals (vdW) semiconductors, such as
transition metal dichalcogenides ($e.g.$ MoS2, WSe2), have appeared as a
fascinating class of materials for exploring novel excitonic phenomena. Mak
_et al._ (2010); Arora  _et al._ (2015a, b); Molas  _et al._ (2017) Inspired
by these exciting achievements, a new group of emerging vdW semiconductors,
$i.e.$ group-IV monochalcogenides MX (where M = Ge, Sn or Pb and X = S, Se or
Te), has attracted increasing attention due to their anisotropic optical
properties. They originate from a low-symmetry orthorhombic crystal structure,
analogous to black phosphorous (BP), which is widely investigated. Ling  _et
al._ (2016a); Ribeiro  _et al._ (2015); Molas  _et al._ (2021) Moreover, the
family of MX materials exhibits high carrier mobility, larger for monolayers
as compared to bulk, Li  _et al._ (2016) which can lead to potential
applications in angle-resolved opto-electronics.

In this work, we investigate anisotropic optical and vibrational properties of
GeS with the aid of the polarization-resolved photoluminescence (PL),
reflectance contrast (RC), and Raman scattering (RS. The low-temperature
($T$=5 K) optical band-gap absorption is governed by a single resonance
related to the neutral exciton, the corresponding emission is dominated by
disorder/impurity- and/or phonon-assisted recombination processes. Both the RC
and PL spectra are found to be linearly polarized along the armchair direction
of the crystal. The low and room ($T$=300 K) temperature RS spectra consist of
six Raman peaks, identified with the help of Density Fuctional Theory (DFT)
calculations. Their polarization properties are studied under four different
excitation energies. We found that the polarization orientations of the
A${}^{2}_{\textrm{g}}$ and A${}^{4}_{\textrm{g}}$ modes under specific
excitation energy can be used to distinguish between armchair and zigzag
crystallographic directions in GeS crystals.

##  II Materials and Methods

###  II.1 Samples

A bulk-like flake of GeS was placed on a Si/(90 nm) SiO2 substrate by
polydimethylsiloxane (PDMS)-based exfoliation of bulk crystals purchased from
HQ Graphene. The flake of interest was initially identified by visual
inspection under an optical microscope then subjected to atomic force
microscopy.

###  II.2 Experimental techniques

The PL spectra measured under laser excitation of $\lambda$= 660 nm (1.88 eV).
The RS measurements were performed using illumination with a series of lasers:
$\lambda$= 488 nm (2.54 eV), $\lambda$= 515 nm (2.41 eV), $\lambda$= 561 nm
(2.21 eV), $\lambda$= 633 nm (1.96 eV). The excitation light in those
experiments was focused by means of a 50x long-working distance objective with
a 0.55 numerical aperture (NA) producing a spot of about 1 $\mu$m diameter.
The signal was collected via the same microscope objective (the backscattering
geometry), sent through a 0.75 m monochromator, and then detected by using a
liquid nitrogen cooled charge-coupled device (CCD) camera. To detect low-
energy RS below 100 cm-1 from the laser line, a set of Bragg filters was
implemented in both excitation and detection paths. In the case of the RC
studies, the only difference in the experimental setup with respect to the one
used for recording the PL and RS signals concerned the excitation source,
which was replaced by a tungsten halogen lamp. The light from the lamp was
coupled to a multimode fiber of a 50 $\mu$m core diameter, and then collimated
and focused on the sample to a spot of about 4 $\mu$m diameter. All
measurements were performed with the samples placed on a cold finger in a
continuous flow cryostat mounted on x–y manual positioners. The excitation
power focused on the sample was kept at 50 $\mu$W during all measurements to
avoid local heating.

The polarization-resolved PL and RC spectra were analyzed by the motorized
half-wave plate and a fixed linear polarizer mounted in the detection path. In
contrast, the polarization-sensitive RS measurements were performed in two co-
(XX) and cross-linear (XY) configurations, which correspond to the parallel
and perpendicular orientation of the excitation and detection polarization
axes, respectively. The analysis of the RS signal was done using a motorized
half-wave plate, mounted on top of the microscope objective, which provides
simultaneous rotation of polarization axis in the XX and XY configurations.

###  II.3 Theoretical calculations

DFT calculations were performed in Vienna Ab initio Simulation Package Kresse
and Furthmüller (1996) with Projector Augmented Wave method Kresse and Joubert
(1999) and Perdew-Burke-Ernzerhof parametrization Perdew  _et al._ (1996) of
general gradients approximation of an exchange-correlation functional. A plane
wave basis energy cutoff of 550 eV and a $\Gamma$-centered Monkhorst-Pack
k-grid of 12$\times$10$\times$4 were found sufficient to converge the lattice
constants up to 0.001 Å. Geometrical parameters were optimized until the
interatomic forces and stress tensor components were lower than 10-5 eV/Å and
0.01 kbar, respectively. The interlayer vdW interactions were taken into
account by the use of Grimme’s D3 correction Grimme  _et al._ (2010). Phonon
dispersion of GeS was calculated within Parliński-Li-Kawazoe method Parlinski
_et al._ (1997), as implemented in Phonopy software Togo and Tanaka (2015).
The 4$\times$4$\times$2 supercells were employed to find the interatomic force
constants within the harmonic approximation. The irreducible representations
of Raman active phonon modes at $\Gamma$ point were determined with the use of
spglib library Togo and Tanaka (2018).

##  III Results

###  III.1 Crystallographic structure

GeS is a layered material, which crystallizes in a distorted orthorhombic
structure (D${}^{1}6_{2h}$), as shown in Figure 1. That form comprising eight
atoms in the primitive unit cell has been proven to be dynamically and
thermally stable at room temperature Wiedemeier and Siemers (1977). The
puckered honeycomb lattice of GeS has an anisotropic crystal structure
characterized by the two orthogonal armchair ($x$) and zigzag ($y$)
directions, denoted in Figure 1(a). The stack of consecutive layers in the
perpendicular direction ($z$) in respect to the $xy$ plane is presented in
Figure 1(b). Note that the crystallographic structure of GeS is analogous to
the BP one.Molas  _et al._ (2021)

Figure 1: (a) Top and side views of the GeS crystal structure for a single
layer. The armchair and zigzag directions are shown in relation to the crystal
orientation and lattice parameters. (b) Geometrical structure of multilayer
GeS.

###  III.2 Optical properties

Figure 2: (a)The low-temperature ($T$=5 K) PL (blue curve) and RC (red curve)
spectra measured on the GeS flake. Polar plots of the integrated intensities
of the X transitions from the (b) PL and (c) RC spectra.  Figure 3: (a) False-
color map of low-temperature PL spectra of GeS as a function of the detection
angle of linear polarization under 1.88 eV laser light excitation. (b)
Temperature evolution of the PL spectra measured on GeS. The spectra are
vertically shifted for clarity and some of them are multiplied by scaling
factors for more clarity.

Figure 2(a) presents the low-temperature ($T$=5 K) PL and RC spectra. The PL
spectrum consists of four emission lines, denoted as X, L1, L2 and L3. In
contrast, the corresponding RC spectrum consists of a single resonance, which
energy coincides with the X emission line. To examine the origin of the X
transition, we carried out the polarization resolved measurements of PL and RC
spectra. The extracted polarization dependencies of the X transitions are
presented in Figures 2(b) and (c). Solid lines in the Figs represent fits of
the evolution of the PL/RC intensity as a function of light polarization. A
dichroic relation of polarized light using Malus law was used:
$I(\theta)$:Hecht (2017)

$I(\theta)=Acos^{2}(\theta-\phi),$ (1)

where $A$ is the amplitude of the emission/absorption transition and $\phi$
represents the phase of polarization dependence. It is seen that both the
emission and absorption signals of the X transitions are linearly polarized
along the same direction (167o and 170o, respectively). The result confirms
directly the same origin of the X feature apparent in both the PL and RC
spectra, as previously reported independently for the emissionHo  _et al._
(2016), photoreflectance Oliva  _et al._ (2020) or absorptionTołłoczko  _et
al._ (2020). Moreover, according to Ref. 20, the X resonance can be related to
the direct transition at the $\Gamma$ point of the Brillouin zone (BZ), which
is linearly polarized along the armchair crystallographic direction.

In order to study the origin of L lines, the polarization-resolved and
temperature evolution of the PL spectra was measured, see Figure 3(a) and (b).
As can be appreciated in the Figure 3(a), the L1, L2 and L3 lines are linearly
polarized along the same armchair direction as the neutral exciton transition.
Moreover, with increasing temperature, the low energy L1, L2 and L3 peaks
quickly disappear from the spectra. At $T$=40 K, only the neutral exciton
contributes to the PL spectrum. The further increase of temperature leads to
the typical redshift and the linewidth broadening of the neutral exciton,
which can be observed up to 120 K. The observed temperature dependence of the
L lines is very similar to the previously reported behavior of so-called
"localized" excitons in monolayers of WS2 and WSe2 exfoliated on Si/SiO2
substrates.Plechinger  _et al._ (2015); Shang  _et al._ (2015); Arora  _et
al._ (2015a); Plechinger  _et al._ (2016); Kłopotowski  _et al._ (2016); Molas
_et al._ (2017) Consequently, we can ascribe tentatively the L1, L2 and L3
peaks to disorder/impurity- and/or phonon-assisted recombination processes.

###  III.3 Vibrational properties

As GeS belongs to the point group D2h, there are 12 Raman-active modes: 4Ag,
2B1g, 2B2g, 4B3g. According to the polarization selection rules for the space
group Pnma (No. 62) of GeS, four Ag and two B1g phonon modes should be
observed in backscattering geometry along the $z$ crystallographic direction.
Their corresponding atomic displacements, denoted by green arrows, are
presented in Figure 4(a). The modes are classified according to their
irreducible representations in the point group D2h, and additionally numbered
due to their increased Raman shift (top index). As can be seen in the Figure
4(a), the Ag modes presents atom movement mostly in the plane defined by the
armchair and out-of-plane directions, while the B1g vibrations take place
along zigzag direction. Figure 4(b) presents the calculated phonon dispersion
with marked Ag and B1g phonon modes active in our experimental conditions. To
verify the theoretical calculations, we measured the RS spectra of GeS at low
($T$=5 K) and room ($T$=300 K) temperatures. Both the RS spectra consist of
six Raman modes: A${}^{1}_{\textrm{g}}$, A${}^{2}_{\textrm{g}}$,
A${}^{3}_{\textrm{g}}$, A${}^{4}_{\textrm{g}}$, B${}^{1}_{\textrm{1g}}$ and
B${}^{2}_{\textrm{1g}}$, which energies are in good agreement with theoretical
calculations shown in Figure 4(b). One notes the characteristic effect of
temperature on the observed peaks. It is seen that both the redshift and the
broadening of Raman peaks scales with the increasing phonon energy, $e.g.$ the
temperature-induced change of the energy and linewidth of the
A${}^{1}_{\textrm{g}}$ mode is much smaller as compared to the
A${}^{3}_{\textrm{g}}$.

Figure 4: (a) Atom displacements (green arrows) for the Raman-active modes.
Axes indicate two view perspectives: armchair, and zigzag directions (b) The
calculated phonon dispersion of GeS. (c) RS spectra measured on GeS at low
($T$=5 K) and room ($T$=300 K) temperatures. The spectra are vertically
shifted for clarity.  Figure 5: Polar plots of the integrated intensities of
the phonon modes measured on GeS at $T$=300 K under 1.96 eV excitation. The
green lines on polar plots are along polarization axes of modes.  Figure 6:
Polar plots of the integrated intensities of the phonon modes
B${}^{2}_{\textbf{1g}}$, A${}^{3}_{\textbf{g}}$, and A${}^{4}_{\textbf{g}}$
measured on GeS at $T$=300 K under 2.54 eV, 2.41 eV, and 2.21 eV excitation.
The green lines on polar plots are along polarization axes of modes.

As varying assignments of Raman peaks in GeS were previously reported,Ho  _et
al._ (2016); Tan  _et al._ (2017); Lam  _et al._ (2018); Ribeiro  _et al._
(2019) we measured the polarization-resolved RS at $T$=300 K under 1.96 eV
excitation. Figure 5 presents polar plots of the integrated intensities as a
function of detection angle for all observed phonon modes in co-linear
configuration (XX). As the corresponding results in cross-linear configuration
(XY) do not add additional value, we focused on XX configuration in our
analysis.Ribeiro  _et al._ (2015) Solid lines represent fits of the modes
intensities as a function of light polarization, I($\theta$), described
by:Ribeiro  _et al._ (2015)

$I(\theta)=(|a|sin^{2}(\theta-\phi)+|c|cos\xi
cos^{2}(\theta-\phi))^{2}+|c|^{2}sin^{2}\xi cos^{4}(\theta-\phi),$ (2)

where $|a|$ and $|c|$ are the amplitudes of the phonon modes, $\phi$
represents the phase of polarization dependence, $\xi$ represents the phase
difference. It is seen that the polarization axes of the
A${}^{1}_{\textrm{g}}$, A${}^{2}_{\textrm{g}}$ and A${}^{4}_{\textrm{g}}$
modes, marked by green lines in the Figure 5, are approximately oriented in
the same direction, $i.e.$ 169o, 168o and 169o, respectively. The direction is
the same as the orientation apparent in PL and RC spectra (167o and 170o,
respectively), which corresponds to the armchair crystallographic direction.
The A${}^{3}_{\textrm{g}}$ mode exhibits different polarization axis of about
112o, which does not match any crystallographic direction. In contrast, the
polarization axes of the B${}^{1}_{\textrm{1g}}$ and B${}^{2}_{\textrm{1g}}$
point to the same direction (123o), which is shifted of about 45o from the
crystallographic directions. In terms of observed symmetries of phonon modes,
the A${}^{2}_{\textrm{g}}$ and A${}^{3}_{\textrm{g}}$ modes display 2-fold
symmetry with an angle period of 180o. In contrast, the
A${}^{1/4}_{\textrm{g}}$ (B${}^{1/2}_{\textrm{1g}}$) mode presents the 4-fold
symmetry with an angle period of 90o and with the different (the same)
intensity of perpendicular arms. Consequently, the polarization axes of the
other phonon modes (except the A${}^{3}_{\textrm{g}}$) can be used to
determine the crystallographic direction, but without the discrimination
between the zigzag and armchair directions. Only the 2-fold symmetry of the
A${}^{2}_{\textrm{g}}$ mode allows to determine the armchair direction of the
crystal.

In order to examine the effect of the excitation energy on the polarization
properties of phonon modes, we performed the polarization-resolved RS
experiments under three more different excitations (2.21 eV, 2.41 eV and 2.54
eV). Due to our experimental limitations, the polarization properties of only
three modes, $i.e.$ B${}^{2}_{\textrm{{1g}}}$, A${}^{3}_{\textrm{g}}$ and
A${}^{4}_{\textrm{g}}$, were analyzed, see Figure 6. Three polarization
characteristics of phonon modes can be distinguished: (i) The
B${}^{1}_{\textrm{g}}$ mode conserves the polarization axis and shape under
excitation with different laser wavelength; (ii) The A${}^{3}_{\textrm{g}}$
mode dramatically changes its polarization axis (compare Figures 5 and 6);
Moreover, its symmetry also changes: from 2-fold symmetry under 1.96 eV
excitation to 4-fold symmetry under other excitations; (iii) In contrast, the
4-fold symmetry of the A${}^{4}_{\textrm{g}}$ mode under 1.96 eV excitation
changes gradually to the 2-fold one under 2.41 eV and 2.54 eV excitations. Its
2-fold symmetry can be useful to determine crystallographic direction of the
crystal. One can conclude that the different excitation energies may affect
significantly the shape between 2- and 4-fold of the A${}^{3}_{\textrm{g}}$
and A${}^{4}_{\textrm{g}}$ modes.

Figure 7: The effect of the excitation energy on the polarization axes of
three phonon modes: B${}^{2}_{\textrm{{1g}}}$, A${}^{3}_{\textrm{g}}$ and
A${}^{4}_{\textrm{g}}$.

The influence of the excitation energies on the polarization axes of three
investigated phonon modes, $i.e.$ B${}^{2}_{\textrm{{1g}}}$,
A${}^{3}_{\textrm{g}}$ and A${}^{4}_{\textrm{g}}$, is summarized in Figure 7.
It can be seen that the polarization axes of the B${}^{2}_{\textrm{{1g}}}$ and
A${}^{4}_{\textrm{g}}$ modes do not change significantly as a function of
excitation energy, whereas the difference between the polarization axes of the
A${}^{3}_{\textrm{g}}$ mode increases of about 34∘ between the 1.93 eV and
2.41 eV excitations. Due to the observed behavior of different peaks, we can
assume that: (i) the polarization properties of the B${}^{2}_{\textrm{{1g}}}$
mode can be used to determine the crystallographic axes in GeS, but without
attribution the zigzag and armchair directions; (ii) the variation of the
polarization axis of the A${}^{3}_{\textrm{g}}$ as a function of the
excitation energy suggests that the electron-phonon coupling may change the
A${}^{3}_{\textrm{g}}$ polarization axis (see Ref. 28 for details); (iii) due
to 2-fold symmetries of the polarization properties (except for 1.96 eV) and
the almost fixed polarization axis of the A${}^{4}_{\textrm{g}}$ mode, its
polarization properties are good to identify the zigzag and armchair
crystallographic directions.

It is important to mention that the observed influence of the excitation
energies on the axes and shape of polarization properties of phonon modes is
very similar to those reported for different anisotropic layered materials,
$e.g.$ BP, ReS2, ReSe2, SnSe1-xSx.Kim  _et al._ (2015); Ling  _et al._
(2016b); Choi  _et al._ (2020); Sriv  _et al._ (2020); Zou  _et al._ (2021)
Moreover, as the effect of thickness on the polarization properties of
different modes in BP has been reported,Kim  _et al._ (2015); Ling  _et al._
(2016b) the outline of further research on thin layers of GeS seems to be
clear.

##  IV Conclusions

We have presented systematic studies of the optical and vibrational properties
of GeS. It has been found that while the low-temperature ($T$=5 K) optical
band-gap absorption is governed by a single resonance related to the neutral
exciton, the corresponding emission is dominated by the disorder/impurity-
and/or phonon-assisted recombination processes. Moreover, both the RC and PL
spectra are found to be linearly polarized along the armchair crystallographic
direction. It is proposed to use the effect to determine crystallographic
direction of GeS. The effect of the excitation energy on the polarization
properties of different phonon modes has been analyzed. It has been shown that
the polarization orientation of the A${}^{2}_{\textrm{g}}$ and
A${}^{4}_{\textrm{g}}$ modes under specific excitation energy can be useful
tools to determine the GeS crystallographic directions: armchair and zigzag.
The strong dependence of the A${}^{3}_{\textrm{g}}$ mode polarization on the
excitation light energy strongly suggests its strong coupling to electronic
excitation of the crystal. We believe that the observations will trigger more
theoretical studies to explain the origin of the electron-phonon interaction.

## Acknowledgements

The work has been supported by the National Science Centre, Poland (Grant No.
2017/27/B/ST3/00205 and 2018/31/B/ST3/02111).

## References

  * Mak  _et al._ (2010) K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Phys. Rev. Lett. 105, 136805 (2010).
  * Arora  _et al._ (2015a) A. Arora, M. Koperski, K. Nogajewski, J. Marcus, C. Faugeras, and M. Potemski, Nanoscale 7, 10421 (2015a).
  * Arora  _et al._ (2015b) A. Arora, K. Nogajewski, M. Molas, M. Koperski, and M. Potemski, Nanoscale 7, 20769 (2015b).
  * Molas  _et al._ (2017) M. R. Molas, K. Nogajewski, A. O. Slobodeniuk, J. Binder, M. Bartos, and M. Potemski, Nanoscale 9, 13128 (2017).
  * Ling  _et al._ (2016a) X. Ling, , S. Huang, E. H. Hasdeo, L. Liang, W. M. Parkin, Y. Tatsumi, A. R. T. Nugraha, A. A. Puretzky, P. M. Das, B. G. Sumpter, D. B. Geohegan, J. Kong, R. Saito, M. Drndic, V. Meunier, and M. S. Dresselhaus, Nano Letters 16, 2260 (2016a).
  * Ribeiro  _et al._ (2015) H. B. Ribeiro, M. A. Pimenta, C. J. S. de Matos, R. L. Moreira, A. S. Rodin, J. D. Zapata, E. A. T. de Souza, and A. H. C. Neto, ACS Nano 9, 4270 (2015).
  * Molas  _et al._ (2021) M. R. Molas, L. Macewicz, A. Wieloszyńska, P. Jakóbczyk, A. Wysmołek, R. Bogdanowicz, and J. B. Jasinski, npj 2D Materials and Applications 5 (2021), 10.1038/s41699-021-00263-8.
  * Li  _et al._ (2016) F. Li, X. Liu, Y. Wang, and Y. Li, Journal of Materials Chemistry C 4, 2155 (2016).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, The Journal of Chemical Physics 132, 154104 (2010).
  * Parlinski  _et al._ (1997) K. Parlinski, Z. Q. Li, and Y. Kawazoe, Phys. Rev. Lett. 78, 4063 (1997).
  * Togo and Tanaka (2015) A. Togo and I. Tanaka, Scripta Materialia 108, 1 (2015).
  * Togo and Tanaka (2018) A. Togo and I. Tanaka, (2018), arXiv:1808.01590 [cond-mat.mtrl-sci] .
  * Wiedemeier and Siemers (1977) H. Wiedemeier and P. A. Siemers, Zeitschrift für anorganische und allgemeine Chemie 431, 299 (1977).
  * Hecht (2017) E. Hecht,  _OPTics, ed. 5 ${}^{t}h$_ (Pearson Education Limited, 2017) p. 347.
  * Ho  _et al._ (2016) C. H. Ho, , and J. X. Li, Material Views 5, 1600814 (2016).
  * Oliva  _et al._ (2020) R. Oliva, T. Woźniak, F. Dybala, A. Tołłoczko, J. Kopaczek, P. Scharoch, and R. Kudrawiec, Phys. Rev. B 101, 235205 (2020).
  * Tołłoczko  _et al._ (2020) A. Tołłoczko, R. Oliva, T. Wożniak, J. Kopaczek, P. Scharoch, and R. Kudrawiec, Material Advances 1, 1886 (2020).
  * Plechinger  _et al._ (2015) G. Plechinger, P. Nagler, J. Kraus, N. Paradiso, C. Strunk, C. Schüller, and T. Korn, physica status solidi (RRL) – Rapid Research Letters 9, 457 (2015).
  * Shang  _et al._ (2015) J. Shang, X. Shen, C. Cong, N. Peimyoo, B. Cao, M. Eginligil, and T. Yu, ACS Nano 9, 647 (2015).
  * Plechinger  _et al._ (2016) G. Plechinger, P. Nagler, A. Arora, A. Granados del Águila, M. V. Ballottin, T. Frank, P. Steinleitner, M. Gmitra, J. Fabian, P. C. M. Christianen, R. Bratschitsch, C. Schüller, and T. Korn, Nano Letters 16, 7899 (2016).
  * Kłopotowski  _et al._ (2016) Ł. Kłopotowski, C. Backes, A. A. Mitioglu, V. Vega-Mayoral, D. Hanlon, J. N. Coleman, V. Y. Ivanov, D. K. Maude, and P. Plochocka, Nanotechnology 27, 425701 (2016).
  * Tan  _et al._ (2017) D. Tan, H. E. Lim, F. Wang, N. B. Mohamed, S. Moun, W. Zhang, Y. Miyauchu, M. Ohfuchi, and K. Matsuda, Nano Research 2, 546 (2017).
  * Lam  _et al._ (2018) D. Lam, K.-S. Chen, J. Kang, X. Liu, and M. C. Hersam, Chemistry of Materials 30, 2245 (2018).
  * Ribeiro  _et al._ (2019) H. B. Ribeiro, S. L. L. M. Ramos, L. Seixas, C. J. S. de Matos, and M. A. Pimenta, Phys. Rev. B 100, 094301 (2019).
  * Zou  _et al._ (2021) B. Zou, Y. Wei, Y. Zhou, D. Ke, X. Zhang, M. Zhang, C.-T. Yip, X. Chen, W. Li, and H. Sun, Nanoscale Horiz. 6, 809 (2021).
  * Kim  _et al._ (2015) J. Kim, J.-U. Lee, J. Lee, H. J. Park, Z. Lee, C. Lee, and H. Cheong, Nanoscale 7, 18708 (2015).
  * Ling  _et al._ (2016b) X. Ling, S. Huang, E. H. Hasdeo, L. Liang, W. M. Parkin, Y. Tatsumi, A. R. T. Nugraha, A. A. Puretzky, P. M. Das, B. G. Sumpter, D. B. Geohegan, J. Kong, R. Saito, M. Drndic, V. Meunier, and M. S. Dresselhaus, Nano Letters 16, 2260 (2016b), pMID: 26963685\.
  * Choi  _et al._ (2020) Y. Choi, K. Kim, S. Y. Lim, J. Kim, J. M. Park, J. H. Kim, Z. Lee, and H. Cheong, Nanoscale Horiz. 5, 308 (2020).
  * Sriv  _et al._ (2020) T. Sriv, T. M. H. Nguyen, Y. Lee, S. Y. Lim, V. Q. Nguyen, K. Kim, S. Cho, and H. Cheong, Scientific Reports 10, 11761 (2020).


f\circ(\mathsf{id}_{A_{j}}\times\upsilon_{j,n})\in\theta_{i}(t)(K^{\prime})\\}\subseteq
L$ $\displaystyle\mbox{definition of }\Gamma$ $\displaystyle\mbox{ iff
}\theta_{i}(t)(K^{\prime})\subseteq L^{\prime}$

The last equivalence in this chain can be checked easily using the definitions
of $K^{\prime}$ and $L^{\prime}$. ∎

######  Proposition 6.

Let $\mathcal{X}=(T_{i},T_{j},\theta_{i},\theta_{j})$ be any choice structure
and $\upsilon=(\upsilon_{i},\upsilon_{j})$ the unique morphism into the
terminal choice structure
$\mathcal{U}=(\Omega_{i},\Omega_{j},\mu_{i},\mu_{j})$. Let $\mathcal{B}_{i}$
be the algebra of $\Omega_{i}$. Then
$\mathcal{B}^{ob}_{\mathcal{X},i}=\\{\upsilon_{i}^{-1}[E]\mid
E\in\mathcal{B}_{i}\\}$. A similar claim holds with $j$ for $i$.

###### Proof.

Let $\mathcal{I}_{i}=\\{\upsilon_{i}^{-1}[E]\mid E\in\mathcal{B}_{i}\\}$ and
$\mathcal{I}_{j}=\\{\upsilon_{j}^{-1}[E]\mid E\in\mathcal{B}_{j}\\}$. Because
taking preimages preserves intersections and complements it follows that
$\mathcal{I}_{i}$ and $\mathcal{I}_{j}$ are algebras. Before we argue that
$\mathcal{B}^{ob}_{\mathcal{X},i}\subseteq\mathcal{I}_{i}$, define

$\mathcal{J}_{j}=\\{(\mathsf{id}_{A_{j}}\times\upsilon_{j})^{-1}[E]\mid
E\mbox{ measurable in }A_{j}\times\Omega_{j}\\}.$

Observe that $A_{j}\times(T_{j},\mathcal{I}_{j})=(A_{j}\times
T_{j},\mathcal{J}_{j})$. One can prove both inclusions of this equality by
considering the generating cylinders in the product spaces
$A_{j}\times(T_{j},\mathcal{I}_{j})$ and $A_{j}\times\Omega_{j}$. Then
consider any act $f\in F(A_{j}\times(T_{j},\mathcal{I}_{j}))=F(A_{j}\times
T_{j},\mathcal{J}_{j})$. We use Lemma 4 from Appendix C.3 to show that there
is some act $f^{\prime}\in F(A_{j}\times\Omega_{j})$ such that

$f=f^{\prime}\circ(\mathsf{id}_{A_{j}}\times\upsilon_{j}).$ (8)

To see that all assumptions of Lemma 4 are satisfied, observe first that by
the definition of $\mathcal{J}_{j}$ the function
$(\mathsf{id}_{A_{j}}\times\upsilon_{j})$ is measurable from $(A_{j}\times
T_{j},\mathcal{J}_{j})$ to $A_{j}\times\Omega_{j}$. Moreover, because $f$ is a
measurable function from $(A_{j}\times T_{j},\mathcal{J}_{j})$ to the discrete
uncertainty space $Z$ we have that $f^{-1}[\\{z\\}]$ is measurable in
$(A_{j}\times T_{j},\mathcal{J}_{j})$ for every $z\in Z$. By the definition of
$\mathcal{J}_{j}$ this means that for every $z\in Z$ there is some measurable
set $E_{z}$ in $A_{j}\times\Omega_{j}$ such that
$f^{-1}[\\{z\\}]=(\mathsf{id}_{A_{j}}\times\upsilon_{j})^{-1}[E_{z}]$.

To prove that $\mathcal{B}^{ob}_{\mathcal{X},i}\subseteq\mathcal{I}_{i}$ we
show that for all finite $K,L\subseteq F(A_{j}\times(T_{j},\mathcal{I}_{j}))$
the set $B^{ob}_{\mathcal{X},i}(K,L)$ is in $\mathcal{I}_{i}$. This, and
together with the analogous statement in which $i$ and $j$ are swapped, shows
that $\mathcal{I}_{i}$ and $\mathcal{I}_{j}$ satisfy the properties from
Definition 3. The claim that
$\mathcal{B}^{ob}_{\mathcal{X},i}\subseteq\mathcal{I}_{i}$ then follows
because $\mathcal{B}^{ob}_{\mathcal{X},i}$ and
$\mathcal{B}^{ob}_{\mathcal{X},j}$ are defined to be the smallest algebras
which satisfy these properties.

Consider any finite $K,L\subseteq F(A_{j}\times(T_{j},\mathcal{I}_{j}))$.
Define the finite sets $K^{\prime},L^{\prime}\subseteq
F(A_{j}\times\Omega_{j})$ such that $K^{\prime}=\\{f^{\prime}\mid f\in K\\}$
and $\\{f^{\prime}\mid f\in L\\}$, where $f^{\prime}$ is defined from $f$ as
in (8). It follows immediately from (8) that

$K=\\{g\circ(\mathsf{id}_{A_{j}}\times\upsilon_{j})\mid g\in
K^{\prime}\\}\quad\mbox{and}\quad
L=\\{g\circ(\mathsf{id}_{A_{j}}\times\upsilon_{j})\mid g\in L^{\prime}\\}.$
(9)

We are then going to show that

$B^{ob}_{\mathcal{X},i}(K,L)=(\Gamma(\mathsf{id}_{A_{j}}\times\upsilon_{j})\circ\theta_{i})^{-1}[B^{K^{\prime}}_{L^{\prime}}],$
(10)

where $B^{K^{\prime}}_{L^{\prime}}=\\{C\in\Gamma(A_{j}\times\Omega_{j})\mid
C(K^{\prime})\subseteq L^{\prime}\\}$ is one of the basic measurable sets that
generates the algebra on $\Gamma(A_{j}\times\Omega_{j})$. Because $\upsilon$
is a morphism of choice structures it follows that
$B^{ob}_{\mathcal{X},i}(K,L)=(\mu_{i}\circ\upsilon_{i})^{-1}[B^{K^{\prime}}_{L^{\prime}}]=\upsilon_{i}^{-1}[\mu_{i}^{-1}[B^{K^{\prime}}_{L^{\prime}}]]$.
From this we obtain that
$B^{ob}_{\mathcal{X},i}(K,L)\in\mathcal{I}_{i}=\\{\upsilon_{i}^{-1}[E]\mid
E\in\mathcal{B}_{i}\\}$ because
$\mu_{i}:\Omega_{i}\to\Gamma(A_{j}\times\Omega_{j})$ is measurable and hence
$\mu_{i}^{-1}[B^{K^{\prime}}_{L^{\prime}}]\in\mathcal{B}_{i}$.

To show (10) we unfold the definitions and see that it amounts to the claim
that $\theta_{i}(t)(K)\subseteq L$ iff
$(\Gamma(\mathsf{id}_{A_{j}}\times\upsilon_{j})(\theta_{i}(t)))(K^{\prime})\subseteq
L^{\prime}$ holds for for all $t\in T_{i}$. Using the definition of $\Gamma$
the right side of the equivalence in this claim unfolds to the statement

$\\{g\in K^{\prime}\mid
g\circ\Gamma(\mathsf{id}_{A_{j}}\times\upsilon_{j})\in\theta_{i}(t)(\\{g\circ\Gamma(\mathsf{id}_{A_{j}}\times\upsilon_{j})\mid
g\in K^{\prime}\\})\\}\subseteq L^{\prime}.$

Using the definitions of $K^{\prime}$ and $L^{\prime}$ and the equalities from
(9) one can easily see that this is equivalent to $\theta_{i}(t)(K)\subseteq
L$.

To prove the other inclusion fix a measurable set $E\in\mathcal{B}_{i}$. We
need to show that $\upsilon_{i}^{-1}[E]\in\mathcal{B}^{ob}_{\mathcal{X},i}$.
By the definition of the algebra $\mathcal{B}_{i}$ on the limit $\Omega_{i}$
it follows that there is some $n$ such that $E=\zeta_{i,n}^{-1}[A]$ for some
$A\subseteq\Omega_{i,n}$. By Lemma 7 we know that
$\upsilon_{i,n}^{-1}[A]\in\mathcal{B}^{ob}_{\mathcal{X},i}$. From (7) we have
that $\upsilon_{i,n}=\zeta_{i,n}\circ\upsilon_{i}$ and thus

$\upsilon_{i,n}^{-1}[A]=(\zeta_{i,n}\circ\upsilon_{i})^{-1}[A]=\upsilon_{i}^{-1}[\zeta_{i,n}^{-1}[A]]=\upsilon_{i}^{-1}[E].$

With this we can conclude that
$\upsilon_{i}^{-1}[E]\in\mathcal{B}^{ob}_{\mathcal{X},i}$. ∎

Because the identity morphism is the unique morphism from the universal choice
structure to itself it follows as a simple corollary from Proposition 6 that
$\mathcal{B}^{ob}_{\mathcal{U},i}$ is equal to the algebra $\mathcal{B}_{i}$
of $\Omega_{i}$. This property has already been observed in other settings and
is called “minimality” in [DT08].

Before we provide the proof of Proposition 1 we would like to mention that the
proposition shows that non-redundancy is the same as the notion of
observability from the theory of coalgebras. In [Jac16] a coalgebra is called
“observable” if the unique morphism from the coalgebra to the terminal
coalgebra is monic. In [Rut00] such coalgebras are called “simple”.

###### Proof of Proposition 1.

For the direction from left to right assume that $\mathcal{X}$ is non-
redundant and consider the unique morphism
$\upsilon=(\upsilon_{i},\upsilon_{j})$ from
$\mathcal{X}=(T_{i},T_{j},\theta_{i},\theta_{j})$ into the universal choice
structure $\mathcal{U}$. To show that $\upsilon_{i}$ is injective consider
$t,t^{\prime}\in T_{i}$ such that $t\neq t^{\prime}$. Because $\mathcal{X}$ is
non-redundant $\mathcal{B}^{ob}_{\mathcal{X},i}$ is separating, and thus there
is some $E\in\mathcal{B}^{ob}_{\mathcal{X},i}$ with $t\in E$ and
$t^{\prime}\notin E$. By Proposition 6 there is then some $E^{\prime}$ that is
measurable in $\Omega_{i}$ such that $E=\upsilon_{i}^{-1}[E^{\prime}]$. But
then it must be the case that $\upsilon_{i}(t)\neq\upsilon_{i}(t^{\prime})$,
because otherwise $\upsilon_{i}(t^{\prime})=\upsilon_{i}(t)\in E^{\prime}$ and
thus $t^{\prime}\in\upsilon_{i}^{-1}[E^{\prime}]=E$. The same reasoning works
for $j$ in place of $i$.

For the other direction assume that the $\upsilon_{i}:T_{i}\to\Omega_{i}$ is
injective. To show that $\mathcal{B}^{ob}_{\mathcal{X},i}$ is separating on
$T_{i}$ consider any $t,t^{\prime}\in T_{i}$ with $t\neq t^{\prime}$. Because
$\upsilon_{i}$ is injective we have
$\upsilon_{i}(t)\neq\upsilon_{i}(t^{\prime})$. The algebra on $\Omega_{i}$,
which is defined as the limit algebra of the cochain
$(\Omega_{i,n},\xi_{i,n})_{n\in\omega}$ for discrete spaces $\Omega_{i,n}$, is
separating because any two distinct infinite sequences have distinct
projections to some finite level $n$, where they can be separated. Thus there
must be a measurable set $E$ in $\Omega_{i}$ such that $t\in E$ and
$t^{\prime}\notin E$. By Proposition 6 the set $\upsilon_{i}^{-1}[E]$ is in
$\mathcal{B}^{ob}_{\mathcal{X},i}$ and obviously this set separates $t$ and
$t^{\prime}$. The same argument can be used with $j$ in place of $i$. ∎

##  Appendix E Proofs for Section 4

###  E.1 Preliminary observations

It is easy to check that $\mathcal{P}$, and hence also $\Pi$ and
${\Pi}^{\heartsuit}$ are functors, that is, they preserve identities and
composition. The construction of the universal preference structure can then
be carried out analogously to the construction for the universal choice
structure given above. In fact it is only required to reprove a variant of
Theorem 5 for the functor $\mathcal{P}$, which is done implicitly in the
proofs of Section 3 from [DT08].

###  E.2 Proposition 2

We now argue that the map $m_{X}:\mathcal{P}X\to\mathcal{C}X$ is injective at
every set $X$. To show this, assume we have two preference relations
$\preccurlyeq$ and $\preccurlyeq^{\prime}$ over $X$ such that
$m_{X}({\preccurlyeq})=m_{X}({\preccurlyeq^{\prime}})$. We need to argue that
$x\preccurlyeq x^{\prime}$ iff $x\preccurlyeq^{\prime}x^{\prime}$ for all
$x,x^{\prime}\in X$. Since the situation is symmetric it suffices to check one
direction. Hence assume that $x\preccurlyeq x^{\prime}$. We want to show that
$x\preccurlyeq^{\prime}x^{\prime}$. It suffices to consider the case where
$x\neq x^{\prime}$ because otherwise $x\preccurlyeq^{\prime}x^{\prime}$
follows because $\preccurlyeq^{\prime}$ is reflexive.

As $x\neq x^{\prime}$ and $x\preccurlyeq x^{\prime}$ it follows that it can
not be the case that $x^{\prime}\preccurlyeq x$, because otherwise there would
be a contradiction with the anti-symmetry of $\preccurlyeq$. As we explain in
Remark 1 this use of anti-symmetry is crucial. As $x\preccurlyeq x^{\prime}$
and not $x^{\prime}\preccurlyeq x$ it follows that $x^{\prime}$ is the only
maximal element of the set $\\{x,x^{\prime}\\}$ in the relation
$\preccurlyeq$. Hence
$m_{X}({\preccurlyeq})(\\{x,x^{\prime}\\})=\\{x^{\prime}\\}$.

By the assumption that $m_{X}({\preccurlyeq})=m_{X}({\preccurlyeq^{\prime}})$
it follows that
$m_{X}({\preccurlyeq^{\prime}})(\\{x,x^{\prime}\\})=\\{x^{\prime}\\}$. But
this is only possible if $x\preccurlyeq^{\prime}x^{\prime}$, which is what we
had to show.

###  E.3 Proposition 3

Proposition 3 states that $\lambda$ is a natural transformation from the
functor $\Pi$ to the functor $\Gamma$. It is easy to check that this reduces
to the claim that $m$ is a natural transformation from $\mathcal{P}$ to
$\mathcal{C}$. This means that we need to show that for every function $f:X\to
Y$ it holds that $m_{X}\circ\mathcal{P}f=\mathcal{C}f\circ m_{Y}$. Note that
here the $Y$ and $X$ are swapped because $\mathcal{P}$ and $\mathcal{C}$ are
contravariant functors.

Fix any function $f:X\to Y$, preference relation
${\preccurlyeq}\in\mathcal{P}Y$ and finite set $K\subseteq X$. We have to show
that

$m_{X}(\mathcal{P}f(\preccurlyeq))(K)=\mathcal{C}f(m_{Y}(\preccurlyeq))(K).$

For the left-to-right inclusion take any $x\in
m_{X}(\mathcal{P}f(\preccurlyeq))(K)$. We need to show that
$x\in\mathcal{C}f(m_{Y}(\preccurlyeq))(K)$. This means we want to show that
$x\in f^{-1}[m_{Y}(\preccurlyeq)(f[K])]\cap K$. Our assumption that $x\in
m_{X}(\mathcal{P}f(\preccurlyeq))(K)$ means that $x$ is a maximal element of
the set $K$ in the order ${\preccurlyeq^{f}}=\mathcal{P}f(\preccurlyeq)$.
Hence $x\in K$ and it remains to show that $x\in
f^{-1}[m_{Y}(\preccurlyeq)(f[K])]$, which means that $f(x)$ is a maximal
element of the set $f[K]$ in the order $\preccurlyeq$. So consider any other
element of $f[K]$, which must be of the form $f(x^{\prime})$ for some
$x^{\prime}\in K$, and assume that $f(x)\preccurlyeq f(x^{\prime})$. We need
to show that then also $f(x^{\prime})\preccurlyeq f(x)$. From equation (1) in
Section 4.1 defining $\preccurlyeq^{f}$ it follows that
$x\preccurlyeq^{f}x^{\prime}$. Then we can use that $x$ is a maximal element
in $\preccurlyeq^{f}$ to conclude that $x^{\prime}\preccurlyeq^{f}x$. Using
(1) again, we then obtain the required $f(x^{\prime})\preccurlyeq f(x)$.

Now consider the right-to-left inclusion. Take any $x\in
f^{-1}[m_{Y}(\preccurlyeq)(f[K])]\cap K$. We need to show that $x\in
m_{X}(\mathcal{P}f(\preccurlyeq))(K)$, which means that $x$ is a maximal
element of the set $K$ in the order
${\preccurlyeq^{f}}=\mathcal{P}f(\preccurlyeq)$. Clearly $x\in K$. To show
that $x$ is a $\preccurlyeq^{f}$-maximal element in $K$ pick any other
$x^{\prime}$ in $K$ with $x\preccurlyeq x^{\prime}$. We need to show that
$x^{\prime}\preccurlyeq x$. From $x\preccurlyeq x^{\prime}$ it follows with
(1) that $f(x)\preccurlyeq f(x^{\prime})$. We now use that $x\in
f^{-1}[m_{Y}(\preccurlyeq)(f[K])]$. From this it follows that $f(x)\in
m_{Y}(\preccurlyeq)(f[K])$. This means that $f(x)$ is maximal in the set
$f[K]$. Because $x^{\prime}\in K$ it holds that also $f(x^{\prime})$ is in the
set $f[K]$. By the maximality of $f(x)$ in $f[K]$ it follows from
$f(x)\preccurlyeq f(x^{\prime})$ that $f(x^{\prime})\preccurlyeq f(x)$. With
(1) we obtain $x^{\prime}\preccurlyeq x$.

###  E.4 Theorem 3

In the proof of Theorem 3 we need a further concept from category theory which
is a generalization of monic morphism. A family of morphism $(f_{j}:Y\to
X_{j})_{j\in J}$ for any index set $J$ is _jointly monic_ if for all further
morphisms $g,h:T\to Y$ it holds that if $f_{j}\circ g=f_{j}\circ h$ for all
$j\in J$ then already $g=h$.

Families of monic morphism are closely related to limits of cochains. Using
the universal property of the limit $X_{\omega}$ with projections $\zeta_{n}$
of a cochain $(X_{n},f_{n})_{n\in\omega}$, it is easy to show that the family
of all projections $(p_{n})_{n\in\omega}$ is jointly monic. Moreover, if we
have another object $T$ with a jointly monic family $(g_{n}:T\to
X_{n})_{n\in\omega}$ such that $g_{n}=f_{n}\circ g_{n+1}$ for all $n$ then the
unique morphism $u:T\to X_{\omega}$ that exists because of the universal
property of the limit $X_{\omega}$ is monic.

To prove Theorem 3, let $\mathcal{U}^{\prime}=(\Omega^{\prime},\mu^{\prime})$
be the universal preference structure, presented as coalgebra for
${\Pi}^{\heartsuit}$. We assume that $\Omega^{\prime}$, $\mu^{\prime}$,
$\Omega^{\prime}_{n}$, $\zeta^{\prime}_{n}$, …are defined in the same way as
the objects $\Omega$, $\mu$, $\Omega_{n}$, $\zeta_{n}$, …are defined in
Section D for the universal choice structure, just using ${\Pi}^{\heartsuit}$
instead of ${\Gamma}^{\heartsuit}$.

Then consider the choice structure
$\lambda({\mathcal{U}^{\prime}})=(\Omega^{\prime},{\lambda}^{\heartsuit}_{\Omega^{\prime}}\circ\mu^{\prime})$
where ${\lambda}^{\heartsuit}$ is a natural transformation from
${\Pi}^{\heartsuit}$ to ${\Gamma}^{\heartsuit}$ that is defined such that it
applies $\lambda$ componentwise. Note that this definition of
$\lambda({\mathcal{U}^{\prime}})$ corresponds to the one given in Section 4.3.
Let $\upsilon$ be the unique morphism from $\lambda({\mathcal{U}^{\prime}})$
to the universal choice structure $\mathcal{U}$ that exists according to
Theorem 2.

The claim of Theorem 3 is that this $\upsilon$ is monic. Because in Theorem 2
$\upsilon$ was obtain using the universal property of the limit $\Omega$ from
the family of approximations
$(\upsilon_{n}:\Omega^{\prime}_{n}\to\Omega_{n})_{n\in\omega}$ it suffices to
show that this family is jointly monic.

Define morphisms
$\delta_{0}=\mathord{!}_{\Omega^{\prime}_{0}}:\Omega^{\prime}_{0}\to\Omega_{0}$
and inductively
$\delta_{n+1}={\Gamma}^{\heartsuit}\delta_{n}\circ{\lambda}^{\heartsuit}_{\Omega^{\prime}_{n}}:\Omega^{\prime}_{n+1}\to\Omega_{n+1}$.
One can show with an induction over $n$ that all these $\delta_{n}$ are monic.
The base case this holds because $\Omega^{\prime}_{0}$ is the terminal object
of $\mathsf{Unc}^{2}$ and in the inductive step we use Proposition 2 and the
fact hat ${\Gamma}^{\heartsuit}$ preserves monics, which we show in Section
C.2. The latter needs that the image of the injective measurable function
$\delta_{n}:\Omega^{\prime}_{n}\to\Omega_{n}$ that is preserved has the
discrete algebra. This is the case because one can show that if $A_{i}$ and
$A_{j}$ are finite then so are all the $\Omega^{\prime}_{n}$.

We then prove by induction on $n$ that

$\upsilon_{n}=\delta_{n}\circ\zeta^{\prime}_{n}.$ (11)

It follows that the $\upsilon_{n}$ are jointly monic because the
$\zeta^{\prime}_{n}$ are projections out of the limit $\Omega^{\prime}$ and
hence jointly monic and the $\delta_{n}$ are all monic.

For the base case, of (11), we have that
$\upsilon_{0}=\delta_{0}\circ\zeta^{\prime}_{0}$ because both morphism map to
the terminal object $\Omega_{0}$ of $\mathsf{Unc}^{2}$.

For the inductive step we use the following computation:

$\displaystyle\upsilon_{n+1}$
$\displaystyle={\Gamma}^{\heartsuit}\upsilon_{n}\circ{\lambda}^{\heartsuit}_{\Omega^{\prime}}\circ\mu^{\prime}$
$\displaystyle\mbox{definition of }\upsilon_{n+1}$
$\displaystyle={\Gamma}^{\heartsuit}(\delta_{n}\circ\zeta^{\prime}_{n})\circ{\lambda}^{\heartsuit}_{\Omega^{\prime}}\circ\mu^{\prime}$
induction hypothesis
$\displaystyle={\Gamma}^{\heartsuit}\delta_{n}\circ{\Gamma}^{\heartsuit}\zeta^{\prime}_{n}\circ{\lambda}^{\heartsuit}_{\Omega^{\prime}}\circ\mu^{\prime}$
$\displaystyle{\Gamma}^{\heartsuit}\mbox{ functor}$
$\displaystyle={\Gamma}^{\heartsuit}\delta_{n}\circ{\lambda}^{\heartsuit}_{\Omega^{\prime}_{n}}\circ{\Pi}^{\heartsuit}\zeta^{\prime}_{n}\circ\mu^{\prime}$
$\displaystyle{\lambda}^{\heartsuit}\mbox{ natural transformation}$
$\displaystyle={\Gamma}^{\heartsuit}\delta_{n}\circ{\lambda}^{\heartsuit}_{\Omega^{\prime}_{n}}\circ\tau^{\prime}_{n+1}\circ\mu^{\prime}$
$\displaystyle\mbox{definition of }\tau^{\prime}_{n+1}$
$\displaystyle={\Gamma}^{\heartsuit}\delta_{n}\circ{\lambda}^{\heartsuit}_{\Omega^{\prime}_{n}}\circ\zeta^{\prime}_{n+1}$
$\displaystyle\mbox{uniqueness property of }\mu^{\prime}$
$\displaystyle=\delta_{n+1}\circ\zeta^{\prime}_{n+1}$
$\displaystyle\mbox{definition of }\delta^{\prime}_{n+1}$


# Probability density functions for photon propagation in a binary (isotropic-
Poisson) statistical mixture with unmatched positives/negatives refractive
indexes

Tiziano Binzoni  Department of Radiology and Medical Informatics, University
Hospital, Geneva, 1211, Switzerland<EMAIL_ADDRESS>Alain Mazzolo
Université Paris-Saclay, CEA, Service d’Études des Réacteurs et de
Mathématiques Appliquées, 91191, Gif-sur-Yvette, France

###### Abstract

The exact homogenized probability density function, for a photon making a step
of length $s$ has been analytically derived for a binary (isotropic-Poisson)
statistical mixture with unmatched refractive indexes. The companions, exact,
homogenized probability density function for a photon to change direction
(“scatter”) with an angle $\vartheta$, and the homogenized albedo, have also
been obtained analytically. These functions also hold even in the case of
negative refractive indexes and allow one to reduce hundreds of MC simulations
of photon propagation in complex binary (isotropic-Poisson) statistical
mixtures, to only one MC simulation, for an equivalent homogeneous medium.
Note, that this is not an approximate approach, but a mathematically
equivalent and exact result. Additionally, some tutorial examples of
homogenized MC simulations are also given.

## Introduction

Monte Carlo (MC) simulations of heterogeneous and disordered media represent
an important subject of study, because the underlying theory may be at the
basis of many physical and natural systems [1, 2, 3, 4, 5, 6, 7, 8, 9, 10].
Due to their complex geometries and statistical properties, these media are
very difficult to tackle, in particular, with analytical techniques. For this
reason, numerical MC simulations represent a “gold standard” allowing to solve
many problems in this domain.

Among these media, binary (isotropic-Poisson) statistical mixtures are of
particular interest, because they allow us to develop algorithms permitting to
treat them as they would be homogeneous media. This equivalence greatly
simplifies MC simulations, since it permits to reduce hundreds of independent
MC simulations (the classical approach), necessary to compute the means of the
wanted parameters, to just one MC simulation (the present approach)[11].

It is in this framework that, in a previous work, we were able to calculate
the homogenized probability density function (pdf) for a photon to make a
single step of length $s$, in a binary (isotropic-Poisson) statistical
mixture[12]. This allows us to treat the medium as homogeneous and the pdf
plays the role of the celebrated Beer-Lambert-Bouguer law in classical MC
simulations. Note that this is not an approximate approach to the MC
simulations, applied to the binary (isotropic-Poisson) statistical mixtures,
but rather a mathematically exact equivalence.

One limitation of the results presented in Ref. [12] is that the binary medium
was treated as if the different bunch of materials composing it had the same
refractive index, which is often not the case in reality. Moreover, the
homogenized phase function, describing the changes in photons direction after
a scattering event, was not assessed at that time. This is a non-trivial
points to be considered, if moreover we want a medium composed with materials
of different refractive indexes. In fact, in this case the photons can also
“scatter” on the surfaces separating the different materials, rendering the
exact determination of the homogenized phase function more complex.

This is why, the above-mentioned reasons prompted us to treat in the present
contribution the following points: 1) Consider a binary (isotropic-Poisson)
statistical mixtures with unmatched (different) refractive indexes; 2) Derive
the homogenized pdf for a photon to make a step of length $s$ (SSF: single
step function); 3) Derive the homogenized pdf for a photon to change direction
(“scatter”) with an angle $\vartheta$ (PF: phase function); 4) For the sake of
completeness and, to the best of our knowledge, the lack of literature on the
topic, the possibility in the previous 3 points to have negative refractive
indexes, is also considered.

To facilitate the reading, we will carry on use the word “homogenized” to
clarify that, e.g., a function or a parameter are the result of the
homogenization procedure, and “classical” when they must be interpreted in the
usual sense. If not explicitly stated, the context should be enough to clarify
this point.

## Theory

### The medium

The characteristics of the medium considered in the present contribution — a
binary (isotropic-Poisson) statistical mixture — have already been presented
in other works [13, 14, 15, 12]. For this reason, we will only provide a brief
summary of the main points necessary to understand the theory presented in the
following sections.

Intuitively, a realization of the random 3D medium, i.e., one among many
possibilities, is obtained by:

  1. 1.

Generating a set of planes with a uniformly distributed random orientation.
This is done is such a way that an arbitrary line intersecting the random
planes is divided in $l_{i}$ random segments (see Fig. 1; dashed arrow A)

Figure 1: Simplified 2D intuitive representation of a realization of a binary
(isotropic-Poisson) statistical mixture. The figure extends to the infinite
plane; for obvious reasons, here only a small region is shown. Continuous
black lines represent the planes. Colored regions represent the tessels of
type a and b. The dash-dotted line, near label C, represents the normal to the
plane. For the dashed arrows A, B and C see text.

of lengths described by the probability density function (pdf).

$\displaystyle p_{L}(l)=\frac{1}{L}e^{-\frac{l}{L}}.$ (1)

Thus, when an arbitrary line intersects the planes, the generated incident
angles, $\vartheta_{i}$, (e.g., Fig. 1; dashed arrow C) are uniformly
distributed over $[0,\frac{\pi}{2}]$, i.e., with pdf

$\displaystyle
p_{\vartheta_{i}}(\vartheta_{i})=\frac{2}{\pi}\mathbb{1}_{\left[0,\frac{\pi}{2}\right]}(\vartheta_{i}),$
(2)

where $\mathbb{1}_{A}(x)$ denotes the indicator function

$\displaystyle\mathbb{1}_{A}(x)=\left\\{\begin{array}[]{lll}1&\leavevmode\nobreak\
\mathrm{if\leavevmode\nobreak\ \leavevmode\nobreak\ }x\in A\\\
0&\leavevmode\nobreak\ \mathrm{if\leavevmode\nobreak\ \leavevmode\nobreak\
}x\notin A\,.\end{array}\right.$ (5)

  2. 2.

To obtain a realization of the binary statistical mixture, the different 3D
regions, delimited by the planes, are randomly “colored” with two colors,
i.e.; with probability $1-P$ to be “red” and a probability P to be “sky blue”
(Fig. 1). This will generate “red” and “sky blue” volumes called, in the
present context, “tessels”. As a consequence, an arbitrary line intersecting
the tessels (see Fig. 1; dashed arrow B) will be divided into $\xi_{i}$ random
segments with pdf

$p_{Tes}(\xi;\sigma)=\sigma
e^{-\sigma\xi};\qquad\sigma\in\\{\sigma_{a},\sigma_{b}\\},$ (6)

where

$\displaystyle\sigma_{a}=\frac{1}{L}P;\qquad\sigma_{b}=\frac{1}{L}(1-P).$ (7)

From now on, parameters referring to “red” and “sky blue” tessels will be
labeled with “a” and “b”, respectively.

  3. 3.

The optical parameters will then be linked to the tessels of type a and b,
i.e.: $\mu_{a_{i}}$, absorption coefficient; $\mu_{s_{i}}$, scattering
coefficient; $\mu_{t_{i}}=\mu_{a_{i}}+\mu_{s_{i}}$ extinction coefficient;
$n_{i}$, refractive index; $g_{i}$, mean cosine (see below); with
$i\in\\{a,b\\}$.

Note that in the present contribution we will always consider
$|n_{a}|>|n_{b}|$. This does not affect the generality of the results.

### Homogenized Single step function $p_{s_{\rm mix}}(.)$

We will now derive the homogenized pdf $p_{s_{\rm mix}}(.)$ for a photon to
make a single step of length $s$, i.e.; the homogenized SSF. In this case, we
have three possibilities for the photon: 1) it reaches a distance $s$ inside
the tessel at the point where it is absorbed; 2) it reaches a distance $s$
inside the tessel at the point where it is scattered; 3) it reaches the tessel
boundary where it is reflected/refracted, because $n_{a}\neq n_{b}$ (note that
this can also be seen as a kind of “scattering” since the photon changes
direction).

In a previous contribution [12] we have demonstrated that:

  * •

The homogenized pdf for a photon to make a step of length $s$, remaining
inside a tessel of type $i\in\\{a,b\\}$ (absorbed or scattered), is
$(\mu_{t_{i}}+\sigma_{i})e^{-(\mu_{t_{i}}+\sigma_{i})s}$.

  * •

The homogenized pdf to make a step that falls on the tessel boundary
(reflected/refracted), of type $i\in\\{a,b\\}$, at distance $s$, is also
$(\mu_{t_{i}}+\sigma_{i})e^{-(\mu_{t_{i}}+\sigma_{i})s}$.

  * •

The probability for a photon to remain inside the tessel $i$, after a single
step, is

$\displaystyle
P_{0}(\mu_{t_{i}},\sigma_{i})=\frac{\mu_{t_{i}}}{\mu_{t_{i}}+\sigma_{i}};\qquad
i\in\\{a,b\\}.$ (8)

These results allow us to express in general (i.e.; simultaneously taking into
account absorption, scattering and reflection/refraction) the pdf for a photon
to make a step of length $s$ in a tessel of type $i$ as

$\displaystyle
p_{s}(s;\mu_{t_{i}},\sigma_{i})=P_{0}(\mu_{t_{i}},\sigma_{i})(\mu_{t_{i}}+\sigma_{i})e^{-(\mu_{t_{i}}+\sigma_{i})s}+[1-P_{0}(\mu_{t_{i}},\sigma_{i})](\mu_{t_{i}}+\sigma_{i})e^{-(\mu_{t_{i}}+\sigma_{i})s}=(\mu_{t_{i}}+\sigma_{i})e^{-(\mu_{t_{i}}+\sigma_{i})s}.$
(9)

Thus, the pdf $p_{s_{\rm mix}}(.)$ can be written as

$\displaystyle p_{s_{\rm
mix}}(s;\mu_{t_{a}},\mu_{t_{b}},\sigma_{a},\sigma_{b},P)=(1-P)p_{s}(s;\mu_{t_{a}},\sigma_{a})+Pp_{s}(s;\mu_{t_{b}},\sigma_{b}),$
(10)

where $1-P$ and $P$ are the probabilities to have a tessel of type a or b,
respectively. Finally, by substituting Eqs. (8) and (9) into Eq. (10), we
obtain the desired homogenized SSF

$\displaystyle p_{s_{\rm
mix}}(s;\mu_{t_{a}},\mu_{t_{b}},L,P)=(1-P)\left(\mu_{t_{a}}+\frac{P}{L}\right)e^{-\left(\mu_{t_{a}}+\frac{P}{L}\right)s}+P\left(\mu_{t_{b}}+\frac{1-P}{L}\right)e^{-\left(\mu_{t_{b}}+\frac{1-P}{L}\right)s}.$
(11)

Given that $p_{s_{\rm mix}}(s;.)$ is a pdf, we have
$\int_{0}^{+\infty}p_{s_{\rm mix}}(s;.)ds=1$.

As already explained in references [16, 17], if the photon starts from a
“fixed” boundary, e.g., from a fixed light source or from the surface that
delimits the medium, the homogenized pdf to make a step of length $s$ can be
obtained as

$\displaystyle p_{s_{\rm fix}}(s;\mu_{t_{a}},\mu_{t_{b}}L,P)$
$\displaystyle=\frac{1-\int_{0}^{s}p_{s_{\rm
mix}}(s^{\prime};\mu_{t_{a}},\mu_{t_{b}},L,P)ds^{\prime}}{\int_{0}^{+\infty}s^{\prime}p_{s_{\rm
mix}}(s^{\prime};\mu_{t_{a}},\mu_{t_{b}},L,P)ds^{\prime}}$
$\displaystyle=\frac{\left(\mu_{t_{a}}+\frac{P}{L}\right)\left(\mu_{t_{b}}+\frac{1-P}{L}\right)\left[(1-P)e^{-\left(\mu_{t_{a}}+\frac{P}{L}\right)s}+Pe^{-\left(\mu_{t_{b}}+\frac{1-P}{L}\right)s}\right]L}{L\left[P\mu_{t_{a}}+(1-P)\mu_{t_{b}}\right]+2P^{2}-2P+1}.$
(12)

Even in this case, we obviously get $\int_{0}^{+\infty}p_{s_{\rm
fix}}(s;.)ds=1$. The validity of $p_{s_{\rm mix}}(.)$ [Eq. (11)] and,
consequently, $p_{s_{\rm fix}}(.)$ [Eq. (12)] has been tested with ad hoc
“gold standard” classical MC simulations (see some examples below).

### Homogenized phase function $p_{\Theta}(\vartheta)$

The derivation of the homogenized PF $p_{\Theta}(\vartheta)$, for positive and
negative refractive indexes, is more complex than that of the SSF. For this
reason, we will first summarize the whole procedure with the help of a
schematic. At the top of Fig. 2,

Figure 2: Schematic summarizing the main steps allowing to derive the PF
$p_{\Theta}(\vartheta)$.

$p_{\Theta}(\vartheta)$ represents the wanted PF, where $\vartheta$
(scattering angle) is the change in the photon direction after a ”scattering”
event. As mentioned in the previous section, we will also consider
reflection/refraction as ”scattering” events, due to the observed induced
changes in the photon’s direction.

Thus, to derive $p_{\Theta}(\vartheta)$, we must answer the following
questions:

  1. 1.

Is the scattering event due to a photon that comes from a tessel a or from b?
This choice is made by exploiting the fact that the photon has probability
$1-P$ to be in a, and $P$ to be in b (Fig. 2; 2nd line).

  2. 2.

Is the photon scattered inside the tessel or is it reflected/refracted on the
tessel boundary (Fig. 2; 3rd line)? This choice is made by knowing the
probability $P_{0}(\mu_{t_{i}},\sigma_{i})$ [Eq. (8)] for a photon to remain
inside the tessel $i\in\\{a,b\\}$ after a single step.

  3. 3.

If there is a reflection, is it total or a partial (Fig. 2; 4th line)? This
choice is made by knowing the probability to have a partial reflection $P_{\rm
partial}$; derived thanks to the critical angle $\vartheta_{c}$. Obviously, if
the photon goes from a tessel of type b to a, this case does not exist
(because when going from a low to a higher refractive index there is no
critical angle).

  4. 4.

Does the photon actually undergo a reflection or a refraction (Fig. 2; 5th
line)? This choice is made by using the probabilities $P_{t_{a\rightarrow
b}}(\vartheta)$, $P_{r_{a\rightarrow b}}(\vartheta)$, $P_{t_{b\rightarrow
a}}(\vartheta)$ and $P_{r_{b\rightarrow a}}(\vartheta)$ that will be derived
in the following sections. Note the dependence of the scattering angle.

Finally, at the end of each tree branch of Fig. 2 (6th line), we can derive
the relatives pdfs, taking into account for the different physical conditions
described in points 1–4, and allowing to build the homogenized
$p_{\Theta}(\vartheta)$, i.e.;

$\displaystyle p_{\Theta}(\vartheta)$ $\displaystyle=\left(1-P\right)\times$
$\displaystyle\Bigl{\\{}P_{0}(\mu_{t_{a}},\sigma_{a})p_{\rm
HG}(\vartheta,g_{a})+\left(1-P_{0}(\mu_{t_{a}},\sigma_{a})\right)\left[P_{\rm
partial}(P_{t_{a\rightarrow
b}}(\vartheta)p_{\vartheta_{0}}(\vartheta)+P_{r_{a\rightarrow
b}}(\vartheta)p_{\vartheta_{1}}(\vartheta))+P_{\rm
total}p_{\vartheta_{2}}(\vartheta)\right]\Bigr{\\}}$
$\displaystyle+\phantom{aaa}P\phantom{aa}\times$
$\displaystyle\Bigl{\\{}P_{0}(\mu_{t_{b}},\sigma_{b})p_{\rm
HG}(\vartheta,g_{b})\Bigr{.}+\left(1-P_{0}(\mu_{t_{b}},\sigma_{b})\right)\left[P_{t_{b\rightarrow
a}}(\vartheta)p_{\vartheta_{3}}(\vartheta)+P_{r_{b\rightarrow
a}}(\vartheta)p_{\vartheta_{4}}(\vartheta)\right]\Bigr{\\}}.$ (13)

The validity of $p_{\Theta}(\vartheta)$ [Eq. (13)] has been tested with ad hoc
“gold standard” classical MC simulations (see some examples below). To avoid a
too heavy formalism, the refractive indexes will never appear as variables in
the functions composing Eq. (13). The context will be sufficient to clarify
this point. The parameters $\mu_{t_{i}}$ and $\sigma_{i}$, $i\in\\{a,b\\}$,
will appear only when necessary.

#### Probability density function $p_{\rm HG}(\vartheta,g_{i})$

For the pdf of a photon to have a scattering angle $\vartheta$, due to a
scattering event inside a tessel $i$, we will use the celebrated Henyey-
Greenstein function[18, 19], applicable in a large number of situations.

$\displaystyle p_{\rm
HG}(\vartheta,g_{i})=\frac{1}{2}\frac{(1-g_{i}^{2})\sin(\vartheta)}{\left[1+g_{i}^{2}-2g_{i}\cos(\vartheta)\right]^{\frac{3}{2}}};\qquad
i\in\\{a,b\\}$ (14)

where $\vartheta\in[0,\pi]$ (Fig. 3),

Figure 3: Example of the change in direction $\vartheta$ of a photon
scattering on a scatterer located inside a tessel.

and $g_{i}$ is the mean cosine. However, other functions can obviously be
used, depending on the particular needs of the reader.

We have seen above that the probability for a photon to occur in this type of
scattering event, described by $p_{HG}(\vartheta,g_{i})$, is
$P_{0}(\mu_{t_{i}},\sigma_{i})$ [Eq. (8)].

#### Probability density function $p_{\vartheta_{0}}(\vartheta)$ (partial
refraction: $a\rightarrow b$, $\vartheta_{i}\leq\vartheta_{c}$)

In the present and in the following sections, we will use the fact that[20]:

  * •

The well-known Snell’s equation

$\displaystyle n_{a}\sin(\vartheta_{i})=n_{b}\sin(\vartheta_{t}),$ (15)

holds, in general, for positive or negative refractive indexes. The practical
consequence is that if one of the two refractive indexes is negative, then the
angle of refraction $\vartheta_{t}$ simply becomes $-\vartheta_{t}$.

  * •

In the case of partial refraction, the critical angle is expressed as

$\displaystyle\vartheta_{c}=\arcsin\left(\left|\frac{n_{b}}{n_{a}}\right|\right).$
(16)

Note that in the present contribution we can have partial refraction only when
we go from a$\rightarrow$b, because we have imposed $|n_{a}|>|n_{b}|$.

  * •

If the photon goes from b$\rightarrow$a, Snell’s equation must be written as

$\displaystyle n_{b}\sin(\vartheta_{i})=n_{a}\sin(\vartheta_{t})$ (17)

(no critical angle in this case).

Thus, in this section we will derive the pdf $p_{\vartheta_{0}}(\vartheta)$ in
the case of partial refraction, i.e.; when the random angle
$\vartheta_{i}\leq\vartheta_{c}$ and the photon goes from a$\rightarrow$b. For
explanatory reasons, we will treat before the case where the refractive
indexes are simultaneously positive or negative. This case is schematically
represented in Fig. 4A.

Figure 4: Four different cases of refraction. A and C are related to the pdf
$p_{\vartheta_{0}}(\vartheta)$. B and D are related to the pdf
$p_{\vartheta_{3}}(\vartheta)$

.

Considered that $\vartheta_{i}$ is uniformly distributed over
$[0,\frac{\pi}{2}]$ [Eq. (2)], the probability to have a refraction/reflection
is

$\displaystyle P_{\rm partial}=\frac{\vartheta_{c}}{\pi/2}.$ (18)

Now, in the case of refraction, the relation that links $\vartheta_{t}$ and
$\vartheta$ is

$\displaystyle\vartheta_{t}=\vartheta_{i}+\vartheta.$ (19)

By substituting Eq. (19) into Eq. (15), and by solving for $\vartheta$, we can
express $\vartheta$ as a function of $\vartheta_{i}$.

$\displaystyle\vartheta=\arcsin\left(\frac{n_{a}}{n_{b}}\sin(\vartheta_{i})\right)-\vartheta_{i}\equiv
g(\vartheta_{i}),$ (20)

and

$\displaystyle\frac{dg(\vartheta)}{d\vartheta}=\frac{n_{a}}{n_{b}}\frac{\cos(\vartheta)}{\sqrt{1-\frac{n_{a}^{2}}{n_{b}^{2}}\sin^{2}(\vartheta)}}-1.$
(21)

By considering two cases, for $\vartheta_{i}=0$ and
$\vartheta_{i}=\vartheta_{c}$, Eq. (20) tells us that $\vartheta$

$\displaystyle\vartheta\in[0,\frac{\pi}{2}-\vartheta_{c}].$ (22)

From Eq. (21) we can see that the equation $dg(\vartheta)/d\vartheta=0$ has no
real solutions, therefore, $dg(\vartheta)/d\vartheta$ is either positive or
negative and $g(\vartheta)$ is an injective function over the interval
$[0,\frac{\pi}{2}-\vartheta_{c}]$. As a result, $p_{\vartheta_{0}}(\vartheta)$
can be derived over this interval using the classical method[21]

$\displaystyle
p_{\vartheta_{0}}(\vartheta)=p_{\vartheta_{c}}\left(g^{-1}(\vartheta)\right)\left|\frac{dg^{-1}(\vartheta)}{d\vartheta}\right|,$
(23)

where

$\displaystyle
p_{\vartheta_{c}}(\vartheta_{i})=\frac{1}{\vartheta_{c}}\mathbb{1}_{[0,\vartheta_{c}]}(\vartheta_{i})$
(24)

is the canonical pdf of $\vartheta_{i}$ over $[0,\vartheta_{c}]$,

$\displaystyle
g^{-1}(\vartheta)=\vartheta_{i}=\arccos\left(\frac{n_{a}-n_{b}\cos(\vartheta)}{\left(n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)\right)^{\frac{1}{2}}}\right)$
(25)

(from Eq. (20)) and

$\displaystyle\frac{dg^{-1}(\vartheta)}{d\vartheta}=\frac{n_{a}n_{b}\cos(\vartheta)-n_{b}^{2}}{n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)}.$
(26)

Finally, by substituting Eqs. (24), (25) and (26) into Eq. (23) we get

$\displaystyle
p_{\vartheta_{0}}(\vartheta)=\frac{1}{\vartheta_{c}}\mathbb{1}_{[0,\frac{\pi}{2}-\vartheta_{c}]}(\vartheta)\left|\frac{n_{a}n_{b}\cos(\vartheta)-n_{b}^{2}}{n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)}\right|.$
(27)

It is maybe useful to highlight again the fact that the case where both
$n_{a}$ and $n_{b}$ are negative is equivalent to the case where both $n_{a}$
and $n_{b}$ are positive.

The last remaining step is to express $p_{\vartheta_{0}}(\vartheta)$ in the
case where one of the refractive indexes is negative. This condition implies
that the relationship existing between the incident and the refracted angle is
(Fig. 4C)

$\displaystyle-\theta_{t}=\theta_{i}-\vartheta.$ (28)

By using Eq. (28) in place of Eq. (19), and by repeating the above procedure,
we obtain

$\displaystyle
p_{\vartheta_{0}}(\vartheta)=\frac{1}{\vartheta_{c}}\mathbb{1}_{[0,\frac{\pi}{2}-\text{sign}(n_{a}n_{b})\vartheta_{c}]}(\vartheta)\left|\frac{n_{a}n_{b}\cos(\vartheta)-n_{b}^{2}}{n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)}\right|.$
(29)

It is interesting to note that Eq. (29) also includes Eq. (27). For this
reason, Eq. (29) will be adopted as a final model for
$p_{\vartheta_{0}}(\vartheta)$.

We have seen that the probability to have a reflection or a refraction is
$P_{\rm partial}$ [Eq. (18)], but Eq. (13) tells us that we also need the
probability to have only a refraction $P_{t_{a\rightarrow b}}(\vartheta)$,
i.e.; the probability to have $p_{\vartheta_{0}}(\vartheta)$.
$P_{t_{a\rightarrow b}}(\vartheta)$ is classically given by the Fresnel’s
formula for non-polarized light[22, 23], but where in the present case the
angles of incidence and refraction must be expressed in terms of $\vartheta$.
So, starting from the Fresnel’s formula and substituting Eqs. (19) and (25),
we can write

$\displaystyle P_{t_{a\rightarrow b}}(\vartheta)$
$\displaystyle=1-\frac{1}{2}\left[\frac{\sin^{2}(\vartheta_{i}-\vartheta_{t})}{\sin^{2}(\vartheta_{i}+\vartheta_{t})}+\frac{\tan^{2}(\vartheta_{i}-\vartheta_{t})}{\tan^{2}(\vartheta_{i}+\vartheta_{t})}\right]$
$\displaystyle=1-\frac{1}{2}\left[\frac{\sin^{2}(-\vartheta)}{\sin^{2}(2\vartheta_{i}+\vartheta)}+\frac{\tan^{2}(-\vartheta)}{\tan^{2}(2\vartheta_{i}+\vartheta)}\right]$
$\displaystyle=1-\frac{1}{2}\left[\csc\left(\vartheta+2\arccos\left(\frac{n_{a}-n_{b}\cos(\vartheta)}{\left(n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)\right)^{\frac{1}{2}}}\right)\right)^{2}\sin(\vartheta)^{2}\right.$
$\displaystyle+\left.\cot\left(\vartheta+2\arccos\left(\frac{n_{a}-n_{b}\cos(\vartheta)}{\left(n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)\right)^{\frac{1}{2}}}\right)\right)^{2}\tan(\vartheta)^{2}\right].$
(30)

Finally, the case for negative refractive indexes ($n_{a}$ or $n_{b}$) can be
obtained by following the same approach, but by using Eq. (28) instead of Eq.
(19), i.e.;

$\displaystyle P_{t_{a\rightarrow b}}(\vartheta)=1-\frac{1}{2}$
$\displaystyle\left\\{\left[\csc\left(\vartheta+2\arccos\left(\frac{n_{a}-n_{b}\cos(\vartheta)}{\left(n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)\right)^{\frac{1}{2}}}\right)\right)\sin(\vartheta)\right]^{2\,\text{sign}(n_{a}n_{b})}\right.$
$\displaystyle+\left.\left[\cot\left(\vartheta+2\arccos\left(\frac{n_{a}-n_{b}\cos(\vartheta)}{\left(n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)\right)^{\frac{1}{2}}}\right)\right)\tan(\vartheta)\right]^{2\,\text{sign}(n_{a}n_{b})}\right\\}.$
(31)

We can see that formally Eq. (31) also includes Eq. (30), and this is why,
from now on, we will only use Eq. (31).

#### Probability density function $p_{\vartheta_{1}}(\vartheta_{1})$ (partial
reflection: $a\rightarrow b$, $\vartheta_{i}\leq\vartheta_{c}$)

As in the previous section, the probability to have
$\theta_{i}\in[0,\vartheta_{c}]$, i.e., a reflection/refraction, obviously
remains $P_{\rm partial}$ [Eq. (18)] .

So, in the case of reflection, Fig. 5

Figure 5: This figure represents the reflection of a photon on a tessel
boundary. This holds in general for positive or negative refractive indexes.
However, if the photon goes from $a\rightarrow b$, then there is a critical
angle $\vartheta_{c}$ for $\theta_{i}$ (not represented), while, if the photon
goes from $b\rightarrow a$, there is no critical angle.

shows that

$\displaystyle\vartheta=\pi-2\vartheta_{i}$ (32)

and, as a consequence, $\vartheta$ is uniformly distributed over the interval

$\displaystyle\vartheta\in[\pi-2\vartheta_{c},\pi].$ (33)

This allows us to immediately express $p_{\vartheta_{1}}(\vartheta)$ as

$\displaystyle
p_{\vartheta_{1}}(\vartheta)=\frac{1}{2\vartheta_{c}}\mathbb{1}_{[\pi-2\vartheta_{c},\pi]}(\vartheta).$
(34)

Similarly to the previous section, we will exploit again Fresnel’s formula to
derive the probability $P_{r}(\vartheta)$ to have Eq. (34), i.e.; the
probability that the photon is reflected in the case
$\vartheta_{i}\leq\vartheta_{c}$. By inserting Eq. (32) into Snell’s law

$\displaystyle\vartheta_{t}=\arcsin\left(\frac{n_{a}}{n_{b}}\sin(\frac{\pi-\vartheta}{2})\right),$
(35)

and by inserting Eqs. (32) and (35) into Fresnel’s formula for reflection, we
obtain

$\displaystyle P_{r_{a\rightarrow b}}(\vartheta)$
$\displaystyle=\frac{1}{2}\left[\frac{\sin^{2}(\vartheta_{i}-\vartheta_{t})}{\sin^{2}(\vartheta_{i}+\vartheta_{t})}+\frac{\tan^{2}(\vartheta_{i}-\vartheta_{t})}{\tan^{2}(\vartheta_{i}+\vartheta_{t})}\right]$
$\displaystyle=\frac{1}{2n_{b}^{2}}\left[n_{b}^{2}+\cos(\vartheta)\left(n_{a}^{2}-n_{b}^{2}+n_{a}^{2}\cos(\vartheta)\right)\right]$
$\displaystyle\times\cot\left(\frac{\vartheta}{2}+\arcsin\left(\left|\frac{n_{a}}{n_{b}}\right|\cos\left(\frac{\vartheta}{2}\right)\right)\right)^{2}\csc\left(\frac{\vartheta}{2}+\arccos\left(\left|\frac{n_{a}}{n_{b}}\right|\cos\left(\frac{\vartheta}{2}\right)\right)\right)^{2}$
(36)

According to Fig. 5 and e.g. Ref.[20] we know that Eq. (36) must hold for
positive or negative refractive indexes. It is sufficient to consider only
their absolute values. This is why $|.|$ appears in the equation.

#### Probability density function $p_{\theta_{2}}(\theta)$ (total reflection:
$a\rightarrow b$, $\vartheta_{i}>\vartheta_{c}$)

In the case where the photon goes from $a\rightarrow b$ and
$\vartheta_{i}>\vartheta_{c}$, we always have total reflection. Thus, using
Eq. (18), the probability to have this condition is

$\displaystyle P_{\rm total}=1-P_{\rm partial}.$ (37)

In this case, the relationship between $\vartheta$ and $\vartheta_{i}$ (Fig.
5) is still described by Eq. (32), however now $\vartheta_{i}$ is uniformly
distributed over the interval

$\displaystyle\theta_{i}\in[\vartheta_{c},\pi/2].$ (38)

As a consequence, $\vartheta$ is uniformly distributed over

$\displaystyle\vartheta\in[0,\pi-2\vartheta_{c}].$ (39)

This allows us to directly express the pdf $p_{\vartheta_{2}}(\vartheta)$ as

$\displaystyle
p_{\vartheta_{2}}(\vartheta)=\frac{1}{\pi-2\vartheta_{c}}\mathbb{1}_{[0,\pi-2\vartheta_{c}]}(\vartheta).$
(40)

#### Probability density function $p_{\vartheta_{3}}(\vartheta)$ (refraction:
$b\rightarrow a$; no critical angle)

Now, let us treat the case where the photon travels from $b\rightarrow a$, and
where the refractive indexes are positive or negative (4B). In this case,
there is no critical angle, and the pdf $p_{\vartheta_{i}}(\vartheta_{i})$ is
expressed as in Eq. (2). Snell’s equation is written as reminded in Eq. (17).
The relationship between the different angles is (Fig. 4B)

$\displaystyle\theta_{t}=\theta_{i}-\theta$ (41)

By applying the same cumbersome procedure reported in Eqs. (20)–(27), we
obtain

$\displaystyle
p_{\vartheta_{3}}(\vartheta)=\frac{1}{\pi}\mathbb{1}_{[0,\frac{\pi}{2}-\vartheta_{c}]}(\vartheta)\left|1-\frac{n_{b}^{2}-n_{a}^{2}}{n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)}\right|.$
(42)

If one of the refractive indexes is negative (Fig. 4D),

$\displaystyle\vartheta_{t}=\vartheta-\vartheta_{i},$ (43)

we obtain

$\displaystyle
p_{\vartheta_{3}}(\vartheta)=\frac{1}{\pi}\mathbb{1}_{[0,\frac{\pi}{2}-\text{sign}(n_{a}n_{b})\vartheta_{c}]}(\vartheta)\left|1-\frac{n_{b}^{2}-n_{a}^{2}}{n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)}\right|.$
(44)

Equation (44) also takes into account Eq. (42).

By following the procedure utilized for Eq. (30), the probability
$P_{t_{b\rightarrow a}}(\vartheta)$ to have Eq. (42), i.e., that the photon is
refracted, is

$\displaystyle P_{t_{b\rightarrow
a}}(\vartheta)=1-\frac{1}{2}\left[\csc\left(\vartheta+2\arccos\left(\frac{n_{b}-n_{a}\cos(\vartheta)}{\left(n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)\right)^{\frac{1}{2}}}\right)\right)^{2}\sin(\vartheta)^{2}\right.$
$\displaystyle+\left.\cot\left(\vartheta+2\arccos\left(\frac{n_{b}-n_{a}\cos(\vartheta)}{\left(n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)\right)^{\frac{1}{2}}}\right)\right)^{2}\tan(\vartheta)^{2}\right],$
(45)

and if one of the refractive indexes is negative (Fig. 4D),

$\displaystyle P_{t_{b\rightarrow a}}(\vartheta)=1-\frac{1}{2}$
$\displaystyle\left\\{\left[\csc\left(\vartheta+2\,\text{sign}(n_{a}n_{b})\arccos\left(\frac{n_{b}-n_{a}\cos(\vartheta)}{\left(n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)\right)^{\frac{1}{2}}}\right)\right)\sin(\vartheta)\right]^{2\,\text{sign}(n_{a}n_{b})}\right.$
$\displaystyle+\left.\left[\cot\left(\vartheta+2\,\text{sign}(n_{a}n_{b})\arccos\left(\frac{n_{b}-n_{a}\cos(\vartheta)}{\left(n_{a}^{2}+n_{b}^{2}-2n_{a}n_{b}\cos(\vartheta)\right)^{\frac{1}{2}}}\right)\right)\tan(\vartheta)\right]^{2\,\text{sign}(n_{a}n_{b})}\right\\}.$
(46)

Equation (46) also includes Eq. (45) .

#### Probability density function $p_{\vartheta_{4}}(\vartheta)$ (reflection:
$b\rightarrow a$; no critical angle)

When the photon travels from $b\rightarrow a$ (no critical angle), the angle
$\vartheta$ (as shown in Fig. 5) is uniformly distributed over

$\displaystyle\vartheta\in[0,\pi]$ (47)

and thus the pdf $p_{\vartheta_{4}}(\vartheta)$ can immediately be written as

$\displaystyle
p_{\vartheta_{4}}(\vartheta)=\frac{1}{\pi}\mathbb{1}_{[0,\pi]}(\vartheta).$
(48)

Once again, by inserting

$\displaystyle\vartheta_{i}=\vartheta_{t}=\frac{\pi-\vartheta}{2}$ (49)

(Fig. 5) in Fresnel’s formula, we obtain

$\displaystyle P_{r_{b\rightarrow a}}(\vartheta)$
$\displaystyle=\frac{1}{2n_{a}^{2}}\left[n_{a}^{2}+\cos(\vartheta)\left(n_{b}^{2}-n_{a}^{2}+n_{b}^{2}\cos(\vartheta)\right)\right]$
$\displaystyle\times\cot\left(\frac{\vartheta}{2}+\arcsin\left(\left|\frac{n_{b}}{n_{a}}\right|\cos\left(\frac{\vartheta}{2}\right)\right)\right)^{2}\csc\left(\frac{\vartheta}{2}+\arccos\left(\left|\frac{n_{b}}{n_{a}}\right|\cos\left(\frac{\vartheta}{2}\right)\right)\right)^{2},$
(50)

that holds for any positive/negative refractive index. It is worth noting that
$P_{r_{b\rightarrow a}}(\vartheta)$ [Eq. (50)] can be derived from
$P_{r_{a\rightarrow b}}(\vartheta)$ [Eq. (36)] by exchanging $n_{a}$ and
$n_{b}$.

#### Homogenized probability density function $p_{\Theta}(\vartheta)$ for the
isotropic medium

Finally, the homogenized pdf $p_{\Theta}(\vartheta)$ for the binary
(isotropic-Poisson) statistical mixture can be found by inserting Eqs. (7) and
(8) into Eq. (13) to obtain

$\displaystyle p_{\Theta}(\vartheta)$
$\displaystyle=\left(1-P\right)\Bigl{\\{}\frac{L\mu_{t_{a}}}{L\mu_{t_{a}}+P}p_{\rm
HG}(\vartheta,g_{a})+\frac{P}{L\mu_{t_{a}}+P}\left[P_{\rm
partial}[P_{t_{a\rightarrow
b}}(\vartheta)p_{\vartheta_{0}}(\vartheta)+P_{r_{a\rightarrow
b}}(\vartheta)p_{\vartheta_{1}}(\vartheta)]+P_{\rm
total}p_{\vartheta_{2}}(\vartheta)\right]\Bigr{\\}}$
$\displaystyle+P\Bigl{\\{}\frac{L\mu_{t_{b}}}{L\mu_{t_{b}}+1-P}p_{\rm
HG}(\vartheta,g_{b})\Bigr{.}+\frac{1-P}{L\mu_{t_{b}}+1-P}\Bigl{.}[P_{t_{b\rightarrow
a}}(\vartheta)p_{\vartheta_{3}}(\vartheta)+P_{r_{b\rightarrow
a}}(\vartheta)p_{\vartheta_{4}}(\vartheta)]\Bigr{\\}}.$ (51)

Given that $p_{\Theta}(\vartheta)$ is a pdf, we have
$\int_{0}^{\pi}p_{\Theta}(\vartheta)d\vartheta=1$, as expected.

### Phase function $p_{\Phi}(\varphi)$ (refraction/reflection)

Equation (51) represents the pdf for the polar angle $\vartheta$. However, in
a MC simulation, we also need the azimuthal angle $\varphi$. Due to the
medium’s construction, which involves using the Henyey-Greenstein pdf and the
isotropic random orientation of the tessel’s boundaries, $\varphi$ is
uniformly distributed over $[0,2\pi]$. For this reason, the pdf for $\varphi$
can be expressed directly as

$\displaystyle
p_{\Phi}(\varphi)=\frac{1}{2\pi}\mathbb{1}_{[0,2\pi]}(\varphi).$ (52)

### Derivation of the albedo for the random medium with mixing statistics

To propagate a photon in a MC simulation, obeying to the pdf
$p_{s}(s;\mu_{t_{i}},\sigma_{i})$ [Eq. (11)], it is necessary to determine
whether the photon is absorbed or scattered at each step. To this aim, we must
calculate the homogenized albedo $\Lambda_{\rm mix}$ for the binary
(isotropic-Poisson) statistical mixture, as described for example in Ref.
[12], i.e.;

$\displaystyle\Lambda_{\rm mix}=\frac{\int_{0}^{\Delta s}p_{s_{\rm
mix}}(s;\mu_{s_{a}},\mu_{s_{b}},L,P)ds}{\int_{0}^{\Delta s}p_{s_{\rm
mix}}(s;\mu_{t_{a}},\mu_{t_{b}},L,P)ds}\approx\frac{(1-P)(\mu_{s_{a}}+\frac{P}{L})+P(\mu_{s_{b}}+\frac{1-P}{L})}{(1-P)(\mu_{t_{a}}+\frac{P}{L})+P(\mu_{t_{b}}+\frac{1-P}{L})},$
(53)

where $\Delta s\ll 1$. Note that if $P=0$ or $P=1$ (only one type of tessel),
we obtain, respectively, the expected classical albedos
$\mu_{s_{a}}/\mu_{t_{a}}$ or $\mu_{s_{b}}/\mu_{t_{b}}$, for homogeneous
media[24].

### MC code for the binary (isotropic-Poisson) statistical mixture

As explained in the Introduction section — and as already explained in the
case of a statistical mixture with matched refractive indexes[12] — the
present approach allows us to treat the binary (isotropic-Poisson) statistical
mixture as a homogeneous medium. For this reason, it is sufficient to replace
the look-up tables for the SSF and the PF, of any classical MC code for photon
propagation, with the look-up tables derived from the proposed Eqs. (11) and
(51), respectively. The classical albedo is replaced by Eq. (53). This, with a
slight difference, i.e., when the photon originates from a fixed boundary
(e.g. the boundary of the homogenized medium) or from a fixed light source,
Eq. (11) must be replaced by Eq. (12)[16, 17]. This is due to the fact that
Eq. (11) is not an exponential law (i.e.; it does not represent the classical
Beer-Lambert-Bouguer law). Below, we will see some explanatory example of the
present homogenized MC approach.

The look-up tables can be derived as usual, e.g., by computing

$\displaystyle\chi_{\rm s_{mix}}=\int_{0}^{s}p_{\rm
s_{mix}}(s^{\prime};\mu_{t_{a}},\mu_{t_{b}},L,P)ds^{\prime};\quad
s\in[0,+\infty[,$ (54) $\displaystyle\chi_{\rm s_{fix}}=\int_{0}^{s}p_{\rm
s_{fix}}(s^{\prime};\mu_{t_{a}},\mu_{t_{b}},L,P)ds^{\prime};\quad
s\in[0,+\infty[$ (55)

and

$\displaystyle\chi_{\Theta}=\int_{0}^{\vartheta}p_{\Theta}(\vartheta^{\prime})d\vartheta^{\prime};\quad\vartheta\in[0,\pi],$
(56)

where $\chi_{\rm s_{mix}}$, $\chi_{\rm s_{fix}}$ and $\chi_{\Theta}$ are
uniformly distributed random variables over $[0,1]$, In practical words, the
look-up tables allow us to obtain random $s$ or $\vartheta$ as a function of
$\chi_{\rm s_{mix}}$, $\chi_{\rm s_{fix}}$ or $\chi_{\Theta}$.

## Explanatory examples

### Example 1

The pdfs $p_{\rm s_{mix}}(s;.)$ and $p_{\Theta}(\vartheta)$ have been tested
by developing a “gold standard” classical MC code, similar to that utilized in
Ref[12], and by taking into account for the unmatched refractive indexes and
different tessel types. Few examples are reported in Figures 6 and 7,

Figure 6: Single step function $p_{\rm s_{mix}}(s;.)$ as a function of
different optical and geometrical parameters. The black lines represent the
relative “gold standard” MC data.  Figure 7: Single step function
$p_{\Theta}(\vartheta)$ as a function of different optical and geometrical
parameters, where $\mu_{a_{b}}=\mu_{s_{b}}=0$. The black lines appearing over
the colored lines represent the relative “gold standard” classical MC data.

where the black lines represent the MC data. A nice agreement can be observed
between the “gold standard” classical MC data and the analytical (homogenized)
models.

### Example 2

Figure 8

Figure 8: Total transmitted light $T_{\rm tot}$ through a semi-infinite slab
of thickness $s_{0}$ as a function of $P$. Vertical bars represent standard
errors.

represents a typical MC simulation that can be performed with the pdfs
obtained in the present contribution, i.e.; the medium — binary (isotropic-
Poisson) statistical mixture — is treated as it would be “homogeneous”. In
this example, the medium geometry is a semi-infinite slab of thickness
$s_{0}$. Tessels of type b are considered to be “holes” of air. The light
source is a pencil beam normally impinging on one of the medium surfaces. The
total transmitted light $T_{\rm tot}$ is detected on the surface opposite to
the light source.

The figure nicely reproduces the expected behavior, i.e., for a small $P$, the
probability 1-P of having vessels of type a with non-null absorption
increases. Thus, fewer photons reach the opposite surface and $T_{\rm tot}$
tends to 0. As $P$ increases, $T_{\rm tot}$ also increases because the
fraction of tessels of type b containing air becomes higher, making it easier
for photons to pass through the slab.

The two cases for $L=30$ mm follow the same behavior, but they go more slowly
to 0 for small $P$. It is interesting to note that, in general, if the $n_{a}$
sign becomes negative while the other parameters remain the same, $T_{\rm
tot}$ decreases.

### Example 3

To demonstrate that the behavior of photons observed in the Example 2 section
is mainly due to the absorption (and scattering) properties of type a tessels,
we repeat the same MC explanatory example, but with
$\mu_{a_{a}}=\mu_{s_{a}}=0$. In this case, only “scattering” effects due to
the presence of the tessels boundaries are taken into account (i.e.; no
photons are absorbed; they are eliminated only by exiting the slab). In this
case, Fig. 9

Figure 9: Same MC simulation as in Example 2 section, but with
$\mu_{t_{a}}=0$. Vertical bars represent standard errors.

nicely shows that the influence of $\mu_{t_{a}}$ observed in 8 disappears. The
curves tend to become symmetric around the center of the $P$ axes.

The minimum $T_{\rm tot}$ values for each curve are found for $P$ values where
the probability to have tessels of type a tends to be equal to the probability
to have tessels of type b. In fact, with this condition, we have the highest
density of a-to-b tessel interfaces inside the slab. This makes it more
difficult for photons to pass through the slab, i.e.; more photons are
reflected on the light source side.

As in Example 2 section, if the $n_{a}$ sign becomes negative, and the other
parameters remain the same, $T_{\rm tot}$ decreases coherently.

## Conclusions

In the present contribution we have derived the homogenized SSF, SF and albedo
for a binary (isotropic-Poisson) statistical mixture with unmatched refractive
indexes. The obtained mathematical tools allow us to reduce hundreds of
complex MC simulations, of photon propagation in binary (isotropic-Poisson)
statistical mixture, to only one simpler MC simulation of an equivalent
homogeneous medium. This homogenized MC simulation is not an approximation,
but is a mathematically equivalent approach, built analytically, giving the
same results as the classical MC.

The tutorial examples presented in the last section already show the potential
interest that can be generated by studying the influence of negative
refractive indexes on photon propagation, or the influence of tessel geometry
on photon scattering behavior. To the best of our knowledge, these findings
have never been reported elsewhere. We hope that these preliminary results
will stimulate future investigations in new domains related to binary
(isotropic-Poisson) statistical mixtures.

## References

  * [1] Barthelemy, P., Bertolotti, J. & Wiersma, D. S.  A lévy flight for light.  _Nature_ 453, 495–498, DOI: https://doi.org/10.1038/nature06948 (2008).
  * [2] Svensson, T. _et al._ Holey random walks: Optics of heterogeneous turbid composites.  _Phys. Rev. E_ 87, 022120, DOI: https://doi.org/10.1103/PhysRevE.87.022120 (2013).
  * [3] Svensson, T. _et al._ Light diffusion in quenched disorder: Role of step correlations.  _Phys. Rev. E_ 89, 022141, DOI: https://doi.org/10.1103/PhysRevE.89.022141 (2014).
  * [4] Davis, A. B. & Marshak, A.  Photon propagation in heterogeneous optical media with spatial correlations: enhanced mean-free-paths and wider-than-exponential free-path distributions.  _Journal of Quantitative Spectroscopy and Radiative Transfer_ 84, 3–34, DOI: https://doi.org/10.1016/S0022-4073(03)00114-6 (2004).
  * [5] Kostinski, A. B. & Shaw, R. A.  Scale-dependent droplet clustering in turbulent clouds.  _Journal of Fluid Mechanics_ 434, 389–398, DOI: https://doi.org/10.1017/S0022112001004001 (2001).
  * [6] Malvagi, F., Byrne, R. N., Pomraning, G. C. & Somerville, R. C. J.  Stochastic radiative transfer in partially cloudy atmosphere.  _Journal of Atmospheric Sciences_ 50, 2146 – 2158, DOI: https://doi.org/10.1175/1520-0469(1993)050<2146:SRTIPC>2.0.CO;2 (1993).
  * [7] Larsen, E. W. & Vasques, R.  A generalized linear boltzmann equation for non-classical particle transport.  _Journal of Quantitative Spectroscopy and Radiative Transfer_ 112, 619–631, DOI: https://doi.org/10.1016/j.jqsrt.2010.07.003 (2011).  2009 International Conference on Mathematics and Computational Methods (M&C 2009).
  * [8] Zuchuat, O., Sanchez, R., Zmijarevic, I. & Malvagi, F.  Transport in renewal statistical media: Benchmarking and comparison with models.  _Journal of Quantitative Spectroscopy and Radiative Transfer_ 51, 689–722, DOI: https://doi.org/10.1016/0022-4073(94)90125-2 (1994).
  * [9] Haran, O., Shvarts, D. & Thieberger, R.  Transport in two-dimensional scattering stochastic media: Simulations and models.  _Phys. Rev. E_ 61, 6183–6189, DOI: https://doi.org/10.1103/PhysRevE.61.6183 (2000).
  * [10] Mercadier, N., Guerin, W., Chevrollier, M. & Kaiser, R.  Lévy flights of photons in hot atomic vapours.  _Nat. Phys._ 5, 602–605, DOI: https://doi.org/10.1038/NPHYS1286 (2009).
  * [11] Pomraning, G. C.  _Linear kinetic theory and particle transport in stochastic mixtures_ (World Scientific, 1991).
  * [12] Binzoni, T. & Mazzolo, A.  Probability density function for random photon steps in a binary (isotropic-poisson) statistical mixture.  _Scientific Reports_ 13, 9887, DOI: https://doi.org/10.1038/s41598-023-36919-2 (2023).
  * [13] Larmier, C., Dumonteil, E., Malvagi, F., Mazzolo, A. & Zoia, A.  Finite-size effects and percolation properties of poisson geometries.  _Phys. Rev. E_ 94, 012130, DOI: https://doi.org/10.1103/PhysRevE.94.012130 (2016).
  * [14] Larmier, C., Zoia, A., Malvagi, F., Dumonteil, E. & Mazzolo, A.  Poisson-box sampling algorithms for three-dimensional markov binary mixtures.  _Journal of Quantitative Spectroscopy and Radiative Transfer_ 206, 70–82, DOI: https://doi.org/10.1016/j.jqsrt.2017.10.020 (2018).
  * [15] Larmier, C., Zoia, A., Malvagi, F., Dumonteil, E. & Mazzolo, A.  Neutron multiplication in random media: Reactivity and kinetics parameters.  _Annals of Nuclear Energy_ 111, 391–406, DOI: https://doi.org/10.1016/j.anucene.2017.09.006 (2018).
  * [16] d’Eon, E.  A reciprocal formulation of nonexponential radiative transfer. 1: Sketch and motivation.  _Journal of Computational and Theoretical Transport_ 47, 84–115, DOI: https://doi.org/10.1080/23324309.2018.1481433 (2018).
  * [17] Binzoni, T. & Martelli, F.  Monte carlo simulations in anomalous radiative transfer: tutorial.  _J. Opt. Soc. Am. A_ 39, 1053–1060, DOI: https://doi.org/10.1364/JOSAA.454463 (2022).
  * [18] Henyey, L. G. & Greenstein, J. L.  Diffuse radiation in the galaxy.  _Astrophys. J._ 93, 70–83 (1941).
  * [19] Binzoni, T., Leung, T. S., Gandjbakhche, A. H., Rüfenacht, D. & Delpy, D. T.  Comment on ‘the use of the henyey–greenstein phase function in monte carlo simulations in biomedical optics’.  _Physics in Medicine & Biology_ 51, L39–L41, DOI: https://doi.org/10.1088/0031-9155/51/22/L01 (2006).
  * [20] Veselago, V. G.  The electrodynamics of substances with simultaneously negative values of $\epsilon$ and $\mu$.  _Soviet Physics Uspekhi_ 10, 509, DOI: https://doi.org/10.1070/PU1968v010n04ABEH003699 (1968).
  * [21] Hsu, H.  _Schaum’s Outline of Probability, Random Variables, and Random Processes, Fourth Edition_ (McGraw-Hill Education, 2020).
  * [22] Born, M. & Wolf, E.  _Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light_ (Cambridge University Press, 1986), 6 edn.
  * [23] Hecht, E.  _Optics_ (Pearson, 2012).
  * [24] Martelli, F., Binzoni, T., Bianco, S. D., Liemert, A. & Kienle, A.  _Light Propagation through Biological Tissue and other Diffusive Media: Theory, Solutions, and Validation_ (SPIE PRESS, Bellingham, Washington USA, 2022).

## Author contributions statement

TB and AM contributed equally to this work. All authors reviewed the
manuscript.

## Data availability

The datasets used and/or analysed during the current study available from the
corresponding author on reasonable request.

## Additional information

Competing financial interests: The authors declare no competing financial
interests.