# DEEP2: Deep Learning Powered De-scattering with Excitation Patterning

Navodini Wijethilake  Center for Advanced Imaging, Faculty of Arts and
Sciences, Harvard University, Cambridge, USA  Department of Electronic and
Telecommunication Engineering, University of Moratuwa, Sri Lanka  Mithunjha
Anandakumar  Center for Advanced Imaging, Faculty of Arts and Sciences,
Harvard University, Cambridge, USA  Department of Electronic and
Telecommunication Engineering, University of Moratuwa, Sri Lanka  Cheng Zheng
Department of Mechanical Engineering, Massachusetts Institute of Technology,
77 Massachusetts Ave., Cambridge, MA 02139, USA.  Laser Biomedical Research
Center, Massachusetts Institute of Technology, 77 Massachusetts Ave.,
Cambridge, MA 02139, USA.  Peter T. C. So  Department of Mechanical
Engineering, Massachusetts Institute of Technology, 77 Massachusetts Ave.,
Cambridge, MA 02139, USA.  Laser Biomedical Research Center, Massachusetts
Institute of Technology, 77 Massachusetts Ave., Cambridge, MA 02139, USA.
Department of Biological Engineering, Massachusetts Institute of Technology,
77 Massachusetts Ave., Cambridge, MA 02139, USA.  Murat Yildirim  Laser
Biomedical Research Center, Massachusetts Institute of Technology, 77
Massachusetts Ave., Cambridge, MA 02139, USA.  Picower Institute for Learning
and Memory, Massachusetts Institute of Technology, 77 Massachusetts Ave.,
Cambridge, MA 02139, USA.  Department of Neuroscience, Cleveland Clinic Lerner
Research Institute, Cleveland, OH 44195, USA.  Dushan N. Wadduwage  Center for
Advanced Imaging, Faculty of Arts and Sciences, Harvard University, Cambridge,
USA<EMAIL_ADDRESS>
###### Abstract

Limited throughput is a key challenge in in-vivo deep-tissue imaging using
nonlinear optical microscopy. Point scanning multiphoton microscopy, the
current gold standard, is slow especially compared to the wide-field imaging
modalities used for optically cleared or thin specimens. We recently
introduced “De-scattering with Excitation Patterning or DEEP”, as a widefield
alternative to point-scanning geometries. Using patterned multiphoton
excitation, DEEP encodes spatial information inside tissue before scattering.
However, to de-scatter at typical depths, hundreds of such patterned
excitations are needed. In this work, we present DEEP2, a deep learning based
model, that can de-scatter images from just tens of patterned excitations
instead of hundreds. Consequently, we improve DEEP’s throughput by almost an
order of magnitude. We demonstrate our method in multiple numerical and
physical experiments including in-vivo cortical vasculature imaging up to four
scattering lengths deep, in alive mice.

## Introduction

Imaging biological structures such as neurons or blood vessels deep inside
scattering tissue is an important yet challenging microscopy problem. For such
in-vivo deep tissue experiments, today, the gold standard is the point-
scanning two- (or three-) photon microscope (PSTPM) [1, 2, 3]. Especially,
PSTPM occupies a unique position in neuroscience, due to its high spatial
resolution, low phototoxicity, and ability to penetrate through deep tissue.
In operation, PSTPM focuses femtosecond laser light at long wavelengths
through scattering tissue to excite fluorescent molecules. Due to long
wavelengths, the excitation light sees little scattering. The microscope then
collects emission fluorescence through the same scattering tissue onto a
detector. This emission light however is at a much shorter wavelength, and
hence encounters significant scattering on its way to the detector. For this
reason, multiple points cannot be excited and resolved at the same time. So
imaging is performed sequentially one point at a time. Thus, PSTPM is
inherently a slow imaging technique. When speed is needed, either resolution
or the imaging field of view should be compromised [4].

Despite the limitations due to emission light scattering, wide-field two- and
three-photon microscopes (WFTPM) have been proposed. For instance, temporal
focusing microscopy (TFM) focuses amplified femtosecond laser pulses
temporally, for depth selective WFTPM. TFM provides excellent optical
sectioning and also penetrates well through scattering tissue due to long
multi-photon excitation wavelengths [5, 6]. Yet, the emission light scatters.
Therefore, in wide-field deep-tissue imaging, some photons are mapped onto
incorrect detector pixels, degrading the SNR and spatial resolution. Few
groups, including us, recently proposed to overcome this limitation by
combining TFM with structured light illumination. Since excitation light
penetrates through the scattering medium, structured illuminations can
modulate the imaging field of view before scattering. Thus, the information
about the structured illuminations can be used to de-scatter TFM images
despite degraded detection. Using this principle, Escobet et al. [7] proposed
“TempoRAl Focusing microscopy with single-pIXel detection (TRAFIX)”. They used
coded excitations and single-pixel detection. We proposed ‘De-scattering by
Excitation Patterning (DEEP)’ [8]. Instead of single-pixel detection, we used
wide-field detection. Single-pixel detection completely relies on coded
illuminations for de-scattering. Wide-field detection, however, retains some
spatial information and only partially depends on coded illumination. Thus,
while TRAFIX needs tens of thousands of illumination patterns to reconstruct a
typical WFTPM image, DEEP requires only hundreds of illumination patterns.
Despite such speed improvement, DEEP is still orders of magnitude slower
compared to other depth-resolved wide-field imaging modalities used for
optically cleared or thin specimens. For instance, modalities such as single
plane illuminated light-sheet microscopy (SPIM) and structured-illumination
microscopy (SIM) are either single shot or few-shot acquisition. Therefore, in
this work, we attempt to further increase the throughput of DEEP by utilizing
prior image information within a deep-learning-based computational imaging
framework.

From a computational imaging point of view, the imaging process that describes
the translation of the ideal image $(x)$ to the observed image $(y)$ is called
the forward model $(f)$. During the forward imaging process, the observed
image $(y)$ is degraded by image noise, low-pass filtering, pixel-value
quantization, sub-sampling, and in our case scattering [9]. An inverse model
should be constructed to map the observed $(y)$ to the ideal expected image
$(x)$. In our original work [8], we used an analytical inverse model that
utilized wavelet sparsity priors. However, recently, for such inverse
problems, deep learning has proven impressively capable. [9, 10, 11]. For
example, deep learning has shown outstanding success in image classification
[12, 13], biomedical imaging [14, 15, 16], segmentation [17], prediction [18],
denoising [19] and other linear and non-linear inverse problems [11, 20, 21,
22]. In our previous work we demonstrated that deep learning can reconstruct
fine biological structures such as dendritic spines from scattered TFM images
[23]. Similarly, in this work we present a deep learning-based inverse
reconstruction method for DEEP-TFM.

Our proposed method is termed DEEP2. DEEP2 consists of a learning-powered
inverse model that can reconstruct a de-scattered image from only 32-DEEP
measurements (instead of 256 in our original work). The model architecture is
inspired by the UNet but modified to include a self-attention mechanism on the
expanding path and a terminating reconstruction block. We trained DEEP2 using
simulated data sets. To simulate training data, we first carefully modeled the
forward imaging process of the DEEP microscope. As a part of the forward
model, we also modeled the optical behavior of scattering tissue using Monte
Carlo simulations. We then used the forward model to simulate DEEP images from
PSTPM-like image stacks (simulated or experimentally acquired). The DEEP2,
trained on the simulated data was finally tested on the experimentally
acquired DEEP images of beads and mouse cortical vasculature. Our results
suggest that the proposed method can reconstruct deep tissue images up to four
scattering lengths deep in mouse cortical vasculature, using only tens of
patterned excitations instead of hundreds.

## Results and Discussion

Figure 1: The proposed deep learning based DEEP imaging method. (A) Optical
schematic of the DEEP-TFM microscope: FS laser - Amplified femtosecond laser;
L1 & L2 - Optical relay; DMD - Digital micro-mirror device; L3 - Excitation
tube lens; Dio - Dichroic mirror; Obj - Objective lens; M - Mirror; L4 -
Emission tube lens; Cam - Camera Detector. (B) Schematic of the proposed
physics based forward model and the deep learning based ’DEEP2’ inverse model.
(C) A representative simulated object. (D) Image of the object in ’C’ at the
final image plane of the microscope. (E) The image in ’D’ after detected using
the EMCCD camera. Note the low photon count and resulting noise.

### Deep-learning Powered DEEP Microscopy

In this work, we propose a deep-learning-based inverse model to reconstruct
images from DEEP (De-scattering with Excitation Patterning) measurements.
Fig.1A shows the optical schematic of the DEEP-TFM microscope. Amplified
femtosecond laser pulses are temporally focused for depth-resolved multi-
photon excitation. In the excitation path, a digital micro-mirror device (DMD)
is placed at a conjugate image plane before the excitation tube lens. The DMD
mirrors project a binary patterned excitation onto the sample at the focal
plane. The emission photons are collected by the objective lens and imaged
onto a camera detector. During the experiment, multiple patterns are used to
generate multiple encoded image measurements. Our goal is to learn an inverse
model to reconstruct de-scattered images using the knowledge about the
excitation patterns and measured images (see Fig.1B). To train such a model,
we first generated simulated training data using a physics-based forward image
model of the DEEP-TFM microscope and the scattering process. We discuss the
forward model, the EMCCD detector model, and the scattering model in detail in
the methods section. Excitation patterns are parameters of the said forward
model (see Fig.1B). We first measured these patterns in a calibration
experiment. Next, PSTPM-like ground-truth images were input to the forward
model to simulate DEEP images (see Fig.1C-E). This allowed us to generate
paired DEEP and PSTPM data to train the inverse model. The inverse model was
based on the seminal UNet architecture. We further modified the UNet by adding
an attention mechanism to the expanding path. Details about the inverse model
are presented in the method section.

To test our method, we experimented on three datasets: (1) a mixture of
fluorescent beads, (2) mouse pyramidal neurons with the dendritic arbor, and
(3) mouse cortical vasculature. We synthetically generated PSTPM-like beads
images for training (see methods). Mouse pyramidal neurons and cortical
vasculature were imaged from multiple animals using PSTPM to generate the
training datasets. We then performed physical DEEP experiments on a mixture of
fluorescent beads and cortical vasculature of anesthetized mice (see methods)
and used our inverse model (trained on simulated data) to successfully de-
scatter experimental DEEP measurements.

In the next two sections, we present numerical reconstruction results on
unseen test data sets for each case; and de-scattering results for
experimental DEEP measurements.

Figure 2: DEEP2 Validation results on the Synthetic Fluorescent Beads. The
beads are in the size of 1.5 um and 6 um, with the intensity of the 1.5 um
beads 5X higher than the 6 um beads. (A) & (D) Two synthetic DEEP-TFM
instances generated using the forward model. Only one patterned image, out of
32 patterned images used as the input for the DEEP2, is visualized. (B) & (E)
Corresponding synthetic ground truths for the (A) & (B) instances. (C) & (F)
DEEP2 reconstruction using the 32 patterned DEEP-TFM synthetic images
corresponding to (A) & (D) instances. The intensity along the lines G and H,
shown in blue and yellow on (A) is visualized in (G) & (H). Similarly, the
intensity along the lines I and J are shown in in (I) & (J) plots. (K) shows
the variance of the 32 patterns.

### Numerical Results on Synthetic Test Data

In this section we present numerical results on synthetic test datasets to
validate our trained inverse model. All experiments used only 32 patterned
illuminations which is nearly an order of magnitude less than the current
state of the art [8].

##### Reconstruction Results on Synthetic Beads Data.

Simulated 3D bead stacks were used as the input to the forward model to obtain
the synthetic DEEP-TFM like instances at 5 scattering lengths. Then, the DEEP2
based inverse model was trained on the synthetic instances and was validated
on similarly generated unseen synthetic beads instances. Fig.2 presents two
validated instances, where the average input is shown in the first column, the
ground truth in the second column, and the reconstruction of the DEEP2 in the
third column. In addition, the intensities along the lateral directions (G, H,
I J in Fig.2), were plotted for further clarification. The DEEP2 successfully
reconstructed the small beads with a radius of 1.5 um.

To obtain the best performing DEEP2, we employed several approaches.First, we
compared the vanilla UNet and the UNet with the scSE module, along with the KL
divergence loss, to obtain the architecture-wise best inverse model. We could
clearly observe that, in the validation, SSIM value was $>0.9$ for the
modified UNet and $>0.8$ for the vanilla UNet. Further, we exploited 4
different loss functions on the modified UNet. We clearly demonstrated that
the SSIM values for the MSE, RMSLE, and Smooth L1 losses were numerically low
compared to the KL Divergence loss. Additionally, for the validation
instances, the SSIM, MSE, and PSNR metric distributions were statistically
assessed. For this, Kruskal Wallis nonparametric test was used to examine
whether each metric is similarly distributed across each loss function. Thus,
the KL Divergence loss was significantly different ($p<0.05$) from the other 3
loss functions for all three metrics.

Figure 3: DEEP2 validation results on the mouse pyramidal neurons with the
dendritic arbor for 2 and 6 scattering length depths. (A) Average of the 32
patterned DEEP-TFM images, generated using the forward model for a single
instance at 2 SLS. (B) The ground truth for that instance, a PSTPM image. (C)
The descattered instance using the DEEP2 with the 32 patterned DEEP-TFM images
at 2 SLS. (D)-(F) Another instance at 2 SLS. (E) is a PSTPM image used as the
ground truth, (F) is the DEEP2 reconstruction.

##### Mouse Pyramidal Neuron Reconstructions.

Figure. 3 shows the validation cohort results on the mouse pyramidal neurons
with a dendritic arbor. We have presented two instances corresponding to 2 and
6 scattering length depths and the respective intensities along the lateral
directions on each instance are presented in the last column. Figure. 3 (A)
and (D) are two separate instances, each representing average of 32 patterned
DEEP-TFM images generated using the forward model with the PSTPM mouse
pyramidal neuron images at 2 scattering lengths. Figure. 3 (B) and (E),
represent the PSTPM images, which were used as the ground truth for the
corresponding (A) and (D) instances. Figure. 3 (C) and (F) are the
reconstructions of DEEP2, that de-scattered the patterned 32 DEEP-TFM images
at 2 scattering lengths. The intensity along the yellow line on (A)-(C) and
(D)-(F) are visualized in (M) and (N) graphs respectively. Similarly, Figure.
3 (G) and (J) show an average of 32 patterned DEEP-TFM images, synthesized
using the forward model, for the same structures shown in (A) and (D) at 6
scattering lengths. Figure. 3 (H) and (K) show corresponding ground truth and
(I) and (L) show the reconstruction of DEEP2, at 6 scattering lengths. The
intensity along the yellow line on (G)-(I) and (J)-(L) are represented on the
(O) and (P) graphs respectively. Thus, we could observe at 2 SLS the DEEP2 was
able to reconstruct spines, fine structures along the dendrite.

Figure 4: DEEP2 validation results on the mouse cortical vasculature
structures at 2 and 6 scattering length depths. (A) Average of the 32
patterned DEEP-TFM images, generated using the forward model for a single
instance at 2 SLS. (B) The ground truth for that instance, a three-photon
Fluorescent microscope image. (C) The descattered instance using the DEEP2
with the 32 patterned DEEP-TFM images at 2 SLS. (D)-(F) Another instance at 2
SLS. (E) is a three-photon Fluorescent microscope image used as the ground
truth, (F) is the DEEP2 reconstruction.

##### Mouse Cortical Vasculature Reconstructions.

The best performing model, i.e. the KL divergence loss with the modified UNet,
as the DEEP2, was used to reconstruct mouse cortical vasculature images from
numerical DEEP data. Figure. 4 shows two validation instances at 2 and 4
scattering lengths and intensity profiles along the highlighted lateral
directions are presented in the last column. Figure. 4 (A) and (D) are the
average of 32 patterned DEEP-TFM images, synthesized using the forward model,
for two separate mouse cortical vasculature instances. The corresponding
ground truths, i.e. the three-photon Fluorescent microscope images, are shown
in Figure. 4 (B) and (E). The DEEP2 reconstruction for the same two instances
at 2 scattering lengths are presented in Figure. 4 (C) and (D). The intensity
profiles along the marked yellow lines on (A)-(C) and (D)-(F) are indicated in
(M) and (N) graphs respectively. Similarly, for the same vasculature
instances, the average over the 32 patterned DEEP-TFM instances synthesized
using the forward model at the 6 scattering lengths are shown in Figure. 4 (G)
and (J). The corresponding ground truths are shown in Figure. 4 (H) and (K)
and the DEEP2 reconstructions at 6 scattering lengths are shown in Figure. 4
(I) and (L). The intensity along the lateral directions marked in yellow on
(G)-(I) and (J)-(L) are represented in Figure. 4 (O) and (P).

Figure 5: DEEP2 testing phase: qualitative evaluation on the fluorescent beads
at 5 scattering lengths. (A) Average of the experimentally acquired 32
patterned DEEP-TFM images. (B) DEEP2 reconstruction. (C) Ground truth. (B)-(F)
The yellow colored box on (A)-(C) images are enlarged for close visualization.
Figure 6: DEEP2 testing phase: qualitative evaluation on the mouse vasculature
structures at 2 and 4 scattering lengths. (A) & (E) Average of 32 DEEP-TFM
images acquired at 2 scattering lengths. (B) & (F) DEEP2 reconstruction at 2
scattering lengths. (C) & (G) DEEP-TFM mathematical inverse algorithm
reconstruction with 255 DEEP-TFM images at 2 scattering lengths. (D) & (H)
DEEP-TFM mathematical inverse algorithm reconstruction with 32 DEEP-TFM images
at 2 scattering lengths. (I) & (M) Average of 32 DEEP-TFM images acquired at 4
scattering lengths. (J) & (N) DEEP2 reconstruction at 2 scattering lengths.
DEEP-TFM mathematical inverse algorithm reconstruction (K) & (O) with 255
DEEP-TFM images. (L) & (P) with 32 DEEP-TFM images.

### Results on Experimental Data

In this section, we present experimental validation of DEEP2 by successfully
de-scattering experimentally acquired DEEP measurements using the inverse
model trained on synthetic data. Experiments used only 32 patterned
illuminations and thus are nearly an order of magnitude faster than the
current state of the art [8].

##### DEEP2 Imaging and Reconstruction of a Mixture of Fluorescent Beads.

The DEEP2 inverse model trained on the artificial beads data at 5 scattering
lengths, was tested on experimentally acquired DEEP-TFM images. In accordance
with the training data, similar sized beads at a similar composition were
imaged through a scattering intralipid layer (resulting 5-6 scattering
lengths) using the DEEP-TFM microscope with 32 patterned excitations. The
images were then input to the trained-inverse model to output the
reconstructions. The results are shown in Figure 5. Figure 5 (A) shows the
average of the 32 patterned DEEP-TFM images and (B) shows the DEEP2
reconstruction. Figure 5 (C) shows the ground truth acquired without the
scattering intralipid layer. A region that consists of beads with radius of 1
um, marked as a yellow box on Figure 5 (A), was closely observed to evaluate
the performance. As shown in the Figure 5 (D), the average of 32 patterned
DEEP-TFM images contained faint signs of those structures. But the proposed
DEEP2could reconstruct those fine beads, eliminating the background noise as
shown in Figure 5 (E). A quantitative evaluation was not performed as the
DEEP-TFM field of view and the ground truth FOV weren’t pixel matched due to
practical experimental reasons.

##### DEEP2 Imaging and Reconstruction of in-vivo Mouse Cortical Vasculature.

Next, the DEEP2 model trained on synthetic cortical vasculature data, was
tested on experimentally acquired DEEP-TFM images of mouse cortical
vasculature (see methods section for DEEP-TFM experimental details). Figure. 6
shows several mouse cortical vascular instances at 2 and 4 scattering-length
depths. Figure. 6 (A) and (E) show the average of the 32 patterned DEEP-TFM
images for two separate FOVs at 2 scattering lengths. Figure. 6 (B) and (F)
are the corresponding DEEP2 reconstruction with 32 DEEP-TFM images. Figure. 6
(C) and (G) are the DEEP-TFM mathematical inverse algorithm reconstruction
with 255 DEEP-TFM images[8]. Figure. 6 (D) and (H) are the DEEP-TFM
mathematical inverse algorithm reconstruction with 32 DEEP-TFM images.
Similarly, Figure. 6 (I) and (M) show the average of the 32 patterned DEEP-TFM
images for two separate FOVs at 4 scattering lengths. Figure. 6 (J) and (N)
are the DEEP2 reconstructions with 32 DEEP-TFM images corresponding to (I) and
(M) respectively. Figure. 6 (K) and (O) are the DEEP-TFM mathematical inverse
algorithm reconstruction with 255 DEEP-TFM images at 4 scattering lengths.
Figure. 6 (L) and (P) are the DEEP-TFM mathematical inverse algorithm
reconstruction with 32 DEEP-TFM images at 4 scattering lengths. The DEEP2
reconstruction (with 32 DEEP-TFM images) showed additional details that were
missing in the inverse algorithm reconstructions with 255 DEEP-TFM images
(compare Figure 5 B, F, J, N to Figure 5 C, G, K, O). However, the average of
raw data (see Figure 5 A, E, I, M) suggests that these structures are real and
not false artifacts generated by the network. The inverse Algorithm failed to
reconstruct with 32 DEEP-TFM images (see Figure 5 D, H, L, P).

## Conclusion

In conclusion, we demonstrate the DEEP2 for in-vivo deep tissue imaging with
fine reconstruction that requires a less number of patterned instances,
compared to our original work [8]. Evaluation on the experimental instances of
in-vivo cortical vasculature imaging up to four scattering lengths strengthens
the applicability of DEEP2 for complex biological structures. One of the
limitations of this work can be identified as the high time consumption for
data generation using the forward model with a comprehensive hyperparameter
tuning. Thus, In future work we will focus on implementing a learning based
forward model that can be integrated with the DEEP2 for optimized
reconstruction.

## Methods

### Forward Model

Figure 7: Forward model Figure 8: Forward model.(A1-2) The xy and xz views of
the excitation PSF. (B1-2) The xy and xz views of the emission PSF. (C1) The
scattering point spread function at a two-scattering-length depth. (C2) The
scattering point spread function at a seven-scattering-length depth. (D1) A
simulated beads object. (D2) The simulated DEEP-TFM image of the object in
’D1’ at a seven-scattering-length depth before detection. Note that the
maximum photon count in the image is close to 5 photons. (D3) The simulated
DEEP-TFM image in ’D2’ detected on the simulated EMCCD camera. The scale bars
in ’A2’ and ’B2’ are $5\mu m$. The scale bars in ’A2’, ’B2’, ’C1’, and ’C2’
are $5\mu m$. The scale bar in ’D3’is $20\mu m$.

##### Model Architecture.

To be used as inputs to our deep learning model for training, we generate
DEEP-TFM image stacks synthetically from the PSTPM image stacks using a
forward imaging model (see Figure 7). Let $exPSF(x,y,z)$ be the 3D excitation
point spread function(PSF) and $H_{t}(x,y)$ be the pattern on the DMD placed
at an image plane at the excitation side. Now, $exPSF(x,y)*_{3D}H_{t}$ is the
3D patterned excitation on the sample. Here $*_{3D}$ is the 3D convolution
operation. The result is multiplied element-wise with the object, $X0(x,y,z)$
to get the excited object cube. Then the resulting 3D excited object is
convolved plane-by-plane with the stack of 2D scattering point spread
functions at each depth ($sPSF(x,y,z)$). This operation is represented by
$*_{2D}$. Each $sPSF$ at a certain depth, represents the effective scattered
light distribution from a point source at that depth. We will introduce and
discuss the sPSFs model in the next section in detail. The resulting 3D cube
is then convolved in 3D with the emission PSF ($emPSF(x,y,z)$). These
operations constitute the following equation.

$Y_{t}(x,y,z)=\;\\{\\{(exPSF(x,y,z)\;*_{3D}\;H_{t}(x,y))\;\circ\;X0(x,y,z)\\}\;*_{2D}\;sPSF(x,y,z)\\}\;*_{3D}\;emPSF(x,y,z)\quad\quad(eq1)$

Finally, we get the image on the detector, by selecting the z-plane that
corresponds to the focal plane ($z_{focal}$) from the last output cube as
shown in the Figure 7. Then the noise is added to the resulting image using
the following noise model discussed in the "EMCCD noise model" section.

$\hat{Y}_{t}(x,y)=f_{EM}(\sim Poiss(Y_{t}(x,y,z_{focal})+\sim
Poiss(\bar{Y}_{Dark}))+\sim Normal(0,\sigma_{Read})\quad\quad(eq2)$

Here, $\sim Poiss(\mu)$ denotes the observations drawn from a Poisson
distribution of mean $\mu$; $\sim Norm(\mu,\sigma)$ denotes the observations
drawn from a Normal distribution of mean $\mu$ and standard deviation
$\sigma$; $\bar{Y}_{Dark}$ is the expected value of dark current of the camera
and $f_{EM}(.)$ models the electron multiplying process of the EMCCD camera.
Figure 8 shows the cross sections of the simulated PSFs and a representative
simulated DEEP-TFM images.

For the beads experiments, fully simulated objects were used as $X0$s. For
Neuronal experiments, experimentally acquired PSTPM images were used as $X0$s
after pre-processing them to match the voxel size of the DEEP-TFM experiment.
During training, the PSTPM images were simulated and used as the output (or
the ground truth) of the network. The object $X0$ was convolved with the
excitation PSF to generate the PSTPM image.

Figure 9: Illustration of the light scattering process in a scattering tissue.

##### Scattering point spread function (sPSF).

In this section we discuss the mathematical model of the scattering point
spread function. We define the scattering point spread function as the
effective photon distribution of a certain z-plane below the tissue surface,
due to a point source placed at that z-plane. Using Monte Carlo simulations of
light transport through scattering tissue from a point source, we model how
the down-stream optics see the light distribution from the plane where the
point source is located. Figure 9 Illustrates how the light scattering process
is modeled using Mote Carlo methods. The surface is treated as the $z=0$ axial
plane. A fluorescent light source is placed at $z=-z_{0}$ plane. The lateral
location of this point source is considered as $x=y=0$. A photon is then
launched from the source, placed at $(0,0,-z_{0})$, in the random direction,
$(u_{x},u_{y},u_{z})$. This photon is propagated for s distance (call this a
hop). Then a scattering event occurs, and the direction of the photon is
changed by a deflection angle, $\theta$, and an azimuth angle, $\psi$ (call
this a spin). Then the photon is hopped and spun again and again until it
reaches the surface. Once it reaches the surface, the photon’s final direction
is recorded. If the direction is within the range captured by the numerical
aperture of the objective, the photon is then backpropagated until it reaches
the original launch plane (i.e. $z=-z_{0}$) and its location on the launch
plane is recorded. When the imaging system is focused at this launch plane the
photon appears to have originated from this location. This process is related
a sufficiently large number of times to generate the sPSF at the launch plane.

In the simulation we need to randomly generate the hop distance (s) and spin
angles ($\theta$ and $\psi$) according to the scattering properties of the
tissue. First, the probability distribution of s is related to the scattering
coefficient of the tissue, $\mu_{s}$. By the definition of $\mu_{s}$ [ref],
the probability of transmission of a photon without encountering a scattering
event after path-length $s$ is given by,

$P(S>s)=e^{-\mu_{s}s}\quad\quad(eq3)$

Therefore, can derive the cumulative distribution function (c.d.f) of $s$,
$F(s)$, and probability density function (p.d.f) of $s$, $f(s)$ as,

$F(s)=P(S\leqslant s)=1-e^{-\mu_{s}s}\quad\quad(eq4)$
$f(s)=\frac{d(F(s))}{dt}=\mu_{s}e^{-\mu_{s}s}\quad\quad(eq5)$

In Mote Carlo methods, $f(s)$ is sampled using a computer generated, uniformly
distributed random number $rnd_{1}$ such that,

$rnd_{1}=F(s)\quad\quad(eq6)$

Therefore, we can derive $s$ as a function of $rnd_{1}$ as,

$s=-ln(1-rnd_{1})/\mu_{s}\quad\quad(eq7)$

Second, the deflection angle, $\theta$, is related to the anisotropy, $g$, of
the tissue according to the Henyey-Greenstein scattering function (HG
function) [ref] that mimics the angular dependence of light scattering by
small particles [ref]. According to the HG function, the p.d.f of
$cos(\theta)$ can be written as,

$f(cos(\theta))=\frac{1}{2}\frac{(1-g^{2})}{(1+g^{2}-2g\;cos(\theta))^{(3/2)}}\quad\quad(eq8)$

From the same principles of sampling $f(cos(\theta))$ using a second random
number $rnd_{2}$ it can be shown that,

$cos(\theta)=\frac{1}{2g}(1+g^{2}-(\frac{1-g^{2}}{1-g+2g\;rnd_{2}})^{2})\quad\quad(eq9)$

Last, the azimuth angle, $\psi$, can be picked at random, i.e. $\psi=rnd_{3}$,
where $rnd_{3}$ is a third computer generated normally distributed random
number. For more details we refer interested readers to Jacques and Wang [24].

##### EMCCD noise model.

Next we need to model the noise added to the image as shown in Eq2. We use an
electron multiplication CCD (EMCCD) camera for detection. Let $Y_{Shot}$ be
the number of electrons $\bar{e}$ from the signal detected before the electron
multiplying (EM) process. Note that the expected value of $Y_{Shot}$ (i.e.
$\bar{Y}_{Shot}$) is equal to $Y_{t}(x,y,z_{focal})$ in Eq2. Let $Y_{Dark}$ be
the number of electrons generated by the dark current of the camera. Dark
current is usually listed in the camera specifications in "$\bar{e}/pixel/s$"
and can be used to calculate the expected value of $Y_{Dark}$ (i.e.
$\bar{Y}_{Dark}$) when the exposure time is known. Both $Y_{Shot}$ and
$Y_{Dark}$ are Poisson distributed. Let, $N_{Read}$ be the read noise of the
camera. The standard deviation of the read noise ($\sigma_{Read}$) is listed
in camera specifications in $\bar{e}$. Read noise can be described by a normal
distribution with zero mean and ($\sigma_{Read}$) standard deviation. Let’s
denote the EM process by a statistical function $f_{EM}(.)$. Then we can write
a formulation for the output signal, $Y_{Out}$ ($=\hat{Y}_{t}(x,y)$ in Eq2) ,
in $\bar{e}$ as,

$Y_{Out}=f_{EM}(Y_{Shot}+Y_{Dark})+N_{Read}\quad\quad(eq10)$

We can then write the mean and variation equations for the signals from eq10
as,

$\bar{Y}_{Out}=g_{EM}(\bar{Y}_{Shot}+\bar{Y}_{Dark})\quad\quad(eq11)$
$\sigma_{Out}^{2}=g_{EM}^{2}F^{2}(\sigma_{Shot}^{2}+\sigma_{Dark}^{2})+\sigma_{Read}^{2}\quad\quad(eq12)$

Here $\bar{Y}$ denotes the mean of the random variable $Y$. $g_{EM}$ is the
EM-gain of the camera and it is the average gain added by the EM process. The
range of $g_{EM}$ is listed in camera specs and the exact value is set during
imaging. $F$ quantifies the noise added by the EM process, $f_{EM}(.)$, and is
discussed in the next subsection.

##### Electron Multiplication Process, $f_{EM}(.)$.

The EM process in practice is noisy and $g_{EM}$ is only the average value of
the EM-gain. This noise is quantified by the excess noise factor (ENF), $F$,
defined in Eq12. We can derive an expression for $F^{2}$ from Eq12 as,

$F^{2}=\frac{\sigma_{Out}^{2}-\sigma_{Read}^{2}}{g_{EM}^{2}(\sigma_{Shot}^{2}+\sigma_{Dark}^{2})}\quad\quad(eq13)$

Here $(\sigma_{Shot}^{2}+\sigma_{Dark}^{2})$ is the variation of the signal
that goes into the EM process, and $(\sigma_{Out}^{2}-\sigma_{Read}^{2})$ is
the variation of the signal that comes out of the EM process (if F is 1 there
is no noise added to by the EM process). Robins and Hadwen[25] derived an
expression for $F$ as,

$F^{2}=\frac{1}{g_{EM}}(\frac{2g_{EM}+\alpha-1}{\alpha+1})=2(g_{EM}-1)g_{EM}^{(-(N+1)/N)}+\frac{1}{g_{EM}}\quad\quad(eq14)$

Here, $N$ is the number of EM gain stages. EM process is treated as Bernoulli
process at each EM stage and $\alpha$ is the probability of an electron
multiplication event happening, i.e. the probability of a success event. The
value of α is usually small and is in the order of 1-2% [ref- Nuvu cam data
sheet]. Considering all input electrons to a stage, the gain for each EM stage
can be described by a Binomial distribution with a probability mass function,

$pmf(x_{Ad};x_{in},\alpha)=Pr(X_{Ad}=x_{Ad})={{x_{in}}\choose{x_{Ad}}}\alpha^{x_{Ad}}(1-\alpha)^{(x_{in}-x_{Ad})}\quad\quad(eq15)$

Here $x_{in}$ is the number of input electrons to the gain stage and $X_{Ad}$
is the number of added electrons through the EM gain stage. The output
electrons from the gain therefore equal to, $X_{Ad}+x_{in}$. $N$ such gains
stages are cascaded to get the final EM output signal $f_{EM}(.)$.

However, simulating a cascade of $pmf(.)$ functions of the form Eq15 to add
noise to every pixel in every simulated image is prohibitively slow. Therefore
for each potential input value to $f_{EM}(.)$ we generated a distribution of
$f_{EM}(.)$ output values. Note that this is doable as the inputs to
$f_{EM}(.)$ are integers that come from Poisson distributed variables (i.e.
$Y_{Shot}+Y_{Dark}$). Then during the simulations, we randomly sampled the
output distribution corresponding to the input value to get the output of
$f_{EM}(.)$ in Eq10.

### Inverse Model

Figure 10: Our proposed architecture: UNet with concurrent channel and spatial
attention mechanism

##### Network Architecture.

In this work, we use a deep learning network, called DEEP2, as the inverse
model. Our proposed DEEP2 is inspired by the UNet[26] architecture with some
modifications as illustrated in Figure 10. DEEP2 consists of 8 convolution
blocks, each block containing 2 convolutional layers followed by batch
normalization and a ReLU activation. In the encoder block, output of each
block is sent through a max pooling layer for downsampling. In the decoder
blocks, an additional block, which we call scSE, is integrated. The importance
of this block is discussed in the next section. The skip connections between
the downsample and upsample parts, enable the model to learn the fine grained
details in the biological structures more efficiently. We exploit a separate
block at the end, which consists of a convolutional layer and a ReLU
activation layer, as the final reconstruction step.

##### UNet with Spatial and channel squeeze and excitation.

In any deep learning model, the main purpose of convolution layer is to learn
to capture all local spatial information within each channel and to create
feature maps jointly encoding the spatial and channel dependencies. To improve
joint encoding of spatial and channel information through convolution, it
requires much effort. Thus, Concurrent Spatial and Channel ‘Squeeze and
Excitation’ (scSE) block is introduced[27]. This improves the feature map by
suppressing the weak features and enhancing the important features. Despite,
recalibrating non-significant features towards zero with scSE may increase the
sparsity in deeper layers as well as reduce the parameter learning [28]. We
denote the performance enhancement of the UNet with scSE over the vanilla UNet
in our experiments.

##### Loss Function.

Our learning algorithms focus on solving the optimization problem given by
structural risk minimization,
$\min_{f}\frac{1}{N}\sum_{i=1}^{N}L_{\theta}(I_{f}\hat{Y}_{t},X0)+\lambda
R(I_{f})$, where $\hat{Y}_{t}$ is the set of DEEP-TFM images given to the
Inverse model $I_{f}$ and the $X0$ is the corresponding ground truth
(unscattered) image. The first term is the empirical risk, which is minimized
while in the second term, $R(I_{f})$ i.e. regularization term, represents the
model complexity and $\lambda$ balances both terms. In the past few years, the
attention of the Machine Learning community is drawn by the loss function,
which provides a significant theoretical and practical value to the
optimization task. Therefore, we have exploited 4 loss functions with the
above described architecture as the inverse model to obtain the optimal
learning algorithm. The least squared error given by Equation 1 is the most
frequently used loss function in regression tasks. Let us assume the output of
the inverse model is $\hat{X0}$,

$Loss_{MSE}=\frac{1}{n}\sum_{t=1}^{n}\left(X0-\hat{X0}\right)^{2}$ (1)

This loss function can converge fast as the gradient of the loss decreases
when the loss approach 0. However, a drawback associated with this loss
function is the high sensitivity for outliers. Root Mean Squared Logarithmic
Error (RMSLE) is more robust to the outliers as it nullifies the effect by
considering the $\log$ value.

$Loss_{RMSLE}=\sqrt{\frac{1}{n}\sum_{t=1}^{n}\left(\log\left(X0\right)-\log\left(\hat{X0}\right)\right)^{2}},$
(2)

Therefore, the above loss function is also utilized in our experiments. Next
we also explored the Smooth L1 loss, which is also more robust for the
outliers compared the mean squared error.

$Loss_{SmoothL1}=\left\\{\begin{array}[]{ll}\frac{1}{n}\sum_{t=1}^{n}\frac{1}{2}\left(X0-\hat{X0}\right)^{2}&\text{
if }\frac{1}{n}\sum_{t=1}^{n}\left|X0-\hat{X0}\right|<1\\\
\frac{1}{n}\sum_{t=1}^{n}\left|X0-\hat{X0}\right|-\frac{1}{2},&\text{
otherwise }\end{array}\right.$ (3)

Smooth L1 loss conditionally chooses either the mean squared loss or the
absolute loss, while combining the advantages of both loss functions.
Nonetheless, we explored the usability of KL Divergence as a loss function for
our inverse algorithm.

$Loss_{KLDivergence}=\frac{1}{n}\sum_{t=1}^{n}X0\cdot\log\left(\frac{X0}{\hat{X0}}\right)$
(4)

KL divergence is identified as a loss function, which can handle distributions
more efficiently despite the normality of the approximation error.

#### Training and Evaluating the Inverse model

Learning rate is 0.001 for training the inverse model, while the momentum and
weight decay are 0.9 and 0.0001 respectively. All the experiments are
conducted with the PyTorch [29] learning framework. All the experiments are
done with the Titan RTX available at FAS cluster.

##### Evaluation Metrics.

We mainly used three metrics to evaluate the performance of our model on the
testing data:(1)Structural Similarity Index (SSIM),(2)Mean Squared Error
(MSE),(3)Peak Signal to Noise Ratio (PSNR)

### Training and Validation Datasets

##### Artificial Fluorescent Beads Mixture Dataset.

For testing our framework, we have an experimentally acquired fluorescent
beads instance using the DEEP-TFM. In addition, we require a large number of
instances to optimize the learning based inverse model. Therefore, Artificial
beads data, in the size of the experimentally acquired DEEP-TFM instance, were
generated to use as the input to the forward model. The 3D volumes with beads
in radius 1.5 um and 6 um were generated along with randomly selected
intensities from a normal distribution with mean 1 and standard deviation 0.1.
The 1.5 radius beads intensities were enhanced by 5 times compared to the 6 um
beads. These artificially generated 3D volumes with beads were given as the
input to the forward model to generate the synthetic DEEP-TFM instances. Since
the experimental beads data acquired using DEEP-TFM was obtained at 5
scattering lengths, the same depth had been used for the synthetic bead
dataset creation. Thus, for 5 scattering length, 3455 instances of were
obtained which we split into 4:1 ratio for training and validation. In
addition, 128 artificial beads instances were generated for independent
testing. Each instance consists of 32 patterned DEEP-TFM like images of size
256 pixels by 256 pixels.

##### Mouse Pyramidal Neuron Dataset.

Experimentally obtained images of pyramidal neurons with the dendritic arbor,
using a PSTPM, were also used for synthetic DEEP-TFM image generation using
the forward model. Each neuronal image volume was 607 pixels * 607 pixels *
263 pixels with a spatial resolution of 0.25 um. We used 11 such neuronal
image volumes for creation of synthetic DEEP-TFM neuronal instances. In order
to enhance the intensity of the fine structures such as dendritic spines, the
intensity was thresholded at 20, by assigning 20 for all the pixels with
intensity over 20. The successive images along the z plane of the image volume
were containing similar structures. Thus, to avoid repeating similar training
instances, the sub image volumes that matches the third dimension of the
exPSF, were extracted from the neuronal volume, to use as the input to the
forward model, with 5 image plane apart from each other. Further, to avoid the
image volumes with no significant structures, the sub image volumes with the
mean exceeding the mean of the whole neuronal image volume were considered as
the input of the forward model. The 3895 synthetic DEEP-TFM neuronal instances
were generated for 2 ( 100 um) and 6 ( 300 um) scattering lengths. Each
instance consists of 326 pixels by 326 pixels ground truth
($\hat{X}_{0}(x,y)$) and 32 patterned images ($t$ = 32) of size 326 pixels by
326 pixels ($\hat{Y}_{t}(x,y)$) as the input of the inverse model $I$.

##### Mouse Cortical Vasculature Dataset.

In addition to the experimentally acquired beads instances, we also acquired
the blood vessel instances at increasing depths using the DEEP-TFM. Hence, to
train the network, we used an equivalent vascular structural image volume
acquired using a three-photon Fluorescent microscope [30]. The vascular image
stack consists of 800 z-stacks, each had an increment of 1.5 um and each sized
512 pixels * 512 pixels with a spatial resolution of 0.6um corresponding to a
field of view of 300um.

##### Artificial Vasculature Dataset.

Additionally, artificial blood vessel volumes were employed [31, 32] to
demonstrate the applicability of this framework on vasculature structures
furthermore. The synthetic volumes originally consisted of blood vessels with
radius ranging from 1-7 pixels. The finest blood vessels with low intensities
were removed by applying a threshold at 190. The blood vessel volumes were
originally sized 325 pixels by 304 pixels by 600 pixels and 20 such volumes
were downloaded. Each volume was rescaled by factor 1.09, normalized by
dividing by the maximum intensity of each volume and extracted volumes of 326
pixels by 326 pixels by 200 pixels volumes to redirect as the input to the
Forward model. Similarly, we generated datasets at 2 and 6 scattering lengths
for blood vessels.

## Acknowledgements

This work was supported by the Center for Advanced Imaging at Harvard
University (DNW, NW, MA), 5-P41EB015871-32 (DNW, PS), R21 MH130067 (PS, DNW),
R21 NS105070 (PS), and R00EB027706 (MY). DNW is also supported by the John
Harvard Distinguished Science Fellowship Program within the FAS Division of
Science of Harvard University. We also thank Dr. Josiah R. Boivin for
providing point scanning two-photon datasets of pyramidal neurons.

## References

  * [1] Rocheleau, J. V. & Piston, D. W.  Two-photon excitation microscopy for the study of living cells and tissues.  _Current protocols in cell biology_ 20, 4–11 (2003).
  * [2] Yildirim, M., Sugihara, H., So, P. T. & Sur, M.  Functional imaging of visual cortical layers and subplate in awake mice with optimized three-photon microscopy.  _Nature communications_ 10, 1–12 (2019).
  * [3] Yildirim, M. _et al._ Label-free three-photon imaging of intact human cerebral organoids for tracking early events in brain development and deficits in rett syndrome.  _eLife_ 11, e78079 (2022).
  * [4] Benninger, R. K. & Piston, D. W.  Two-photon excitation microscopy for the study of living cells and tissues.  _Current protocols in cell biology_ 59, 4–11 (2013).
  * [5] Vaziri, A. & Shank, C. V.  Ultrafast widefield optical sectioning microscopy by multifocal temporal focusing.  _Optics express_ 18, 19645–19655 (2010).
  * [6] Rowlands, C. J. _et al._ Wide-field three-photon excitation in biological samples.  _Light: Science & Applications_ 6, e16255–e16255 (2017).
  * [7] Escobet-Montalbán, A. _et al._ Wide-field multiphoton imaging through scattering media without correction.  _Science advances_ 4, eaau1338 (2018).
  * [8] Zheng, C. _et al._ De-scattering with excitation patterning enables rapid wide-field imaging through scattering media.  _Science Advances_ 7, eaay5496 (2021).
  * [9] Belthangady, C. & Royer, L. A.  Applications, promises, and pitfalls of deep learning for fluorescence image reconstruction.  _Nature methods_ 1–11 (2019).
  * [10] Weigert, M. _et al._ Content-aware image restoration: pushing the limits of fluorescence microscopy.  _Nature methods_ 15, 1090–1097 (2018).
  * [11] Jin, K. H., McCann, M. T., Froustey, E. & Unser, M.  Deep convolutional neural network for inverse problems in imaging.  _IEEE Transactions on Image Processing_ 26, 4509–4522 (2017).
  * [12] Krizhevsky, A., Sutskever, I. & Hinton, G. E.  Imagenet classification with deep convolutional neural networks.  In _Advances in neural information processing systems_ , 1097–1105 (2012).
  * [13] Ziletti, A., Kumar, D., Scheffler, M. & Ghiringhelli, L. M.  Insightful classification of crystal structures using deep learning.  _Nature communications_ 9, 1–10 (2018).
  * [14] Wei, Z., Liu, D. & Chen, X.  Dominant-current deep learning scheme for electrical impedance tomography.  _IEEE Transactions on Biomedical Engineering_ 66, 2546–2555 (2019).
  * [15] Liu, S. _et al._ Early diagnosis of alzheimer’s disease with deep learning.  In _2014 IEEE 11th international symposium on biomedical imaging (ISBI)_ , 1015–1018 (IEEE, 2014).
  * [16] Wang, S. _et al._ Accelerating magnetic resonance imaging via deep learning.  In _2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI)_ , 514–517 (IEEE, 2016).
  * [17] Girshick, R., Donahue, J., Darrell, T. & Malik, J.  Rich feature hierarchies for accurate object detection and semantic segmentation (pp. 580–587).  In _CVPR’14: Proceedings of the 2014 IEEE Conference on Computer Vision and Pattern Recognition, IEEE Computer Society_ (2014).
  * [18] Nielsen, A. A. & Voigt, C. A.  Deep learning to predict the lab-of-origin of engineered dna.  _Nature communications_ 9, 1–10 (2018).
  * [19] Eraslan, G., Simon, L. M., Mircea, M., Mueller, N. S. & Theis, F. J.  Single-cell rna-seq denoising using a deep count autoencoder.  _Nature communications_ 10, 1–14 (2019).
  * [20] Eulenberg, P. _et al._ Reconstructing cell cycle and disease progression using deep learning.  _Nature communications_ 8, 1–6 (2017).
  * [21] Wei, Z. & Chen, X.  Physics-inspired convolutional neural network for solving full-wave inverse scattering problems.  _IEEE Transactions on Antennas and Propagation_ 67, 6138–6148 (2019).
  * [22] Zhu, B., Liu, J. Z., Cauley, S. F., Rosen, B. R. & Rosen, M. S.  Image reconstruction by domain-transform manifold learning.  _Nature_ 555, 487–492 (2018).
  * [23] Wei, Z. _et al._ 3d deep learning enables fast imaging of spines through scattering media by temporal focusing microscopy.  _arXiv preprint arXiv:2001.00520_ (2019).
  * [24] Jacques, S. L. & Wang, L.  Monte carlo modeling of light transport in tissues.  In _Optical-thermal response of laser-irradiated tissue_ , 73–100 (Springer, 1995).
  * [25] Robbins, M. S. & Hadwen, B. J.  The noise performance of electron multiplying charge-coupled devices.  _IEEE transactions on electron devices_ 50, 1227–1232 (2003).
  * [26] Ronneberger, O., Fischer, P. & Brox, T.  U-net: Convolutional networks for biomedical image segmentation.  In _International Conference on Medical image computing and computer-assisted intervention_ , 234–241 (Springer, 2015).
  * [27] Roy, A. G., Navab, N. & Wachinger, C.  Concurrent spatial and channel ‘squeeze & excitation’in fully convolutional networks.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , 421–429 (Springer, 2018).
  * [28] Uhrig, J. _et al._ Sparsity invariant cnns.  In _2017 international conference on 3D Vision (3DV)_ , 11–20 (IEEE, 2017).
  * [29] Paszke, A. _et al._ Automatic differentiation in pytorch.  _-_ (2017).
  * [30] Yildirim, M. _et al._ Quantitative third-harmonic generation imaging of mouse visual cortex areas reveals correlations between functional maps and structural substrates.  _Biomedical Optics Express_ 11, 5650–5673 (2020).
  * [31] Schneider, M., Reichold, J., Weber, B., Székely, G. & Hirsch, S.  Tissue metabolism driven arterial tree generation.  _Medical image analysis_ 16, 1397–1414 (2012).
  * [32] Todorov, M. I. _et al._ Machine learning analysis of whole mouse brain vasculature.  _Nature methods_ 17, 442–449 (2020).


# Electromagnetic Energy Density in Hyperbolic Metamaterials

Afshin Moradi1<EMAIL_ADDRESS>Pi-Gang Luan2<EMAIL_ADDRESS>_1
Department of Engineering Physics, Kermanshah University of Technology,
Kermanshah, Iran
2Department of Optics and Photonics, National Central University, Jhongli
District, Taoyuan City 320, Taiwan _

###### Abstract

We derive the energy density associated with an electromagnetic wave passing
through a hyperbolic metamaterial (HMM). Both types of HMMs are studied. By
considering a dispersive and absorbing HMM as an effective uniaxial crystal,
we find that under the influence of an electromagnetic wave, the responses of
type I (type II) HMM in the directions perpendicular to and parallel to the
optical axis are similar to those of Lorentz (Drude) and Drude (Lorentz)
media, respectively. Numerical examples are presented to reveal the general
characteristics of the direction-dependent energy storage capacity of both
types of HMMs.

###### pacs:

41.20.Jb,

##  I Introduction

Metamaterials are artificial media that usually refers to arrays of wires and
split-ring resonators (SRRs) D.R.S4184  and have unusual optical properties
such as negative refractive index J.B.P3966 , subwavelength imaging X.Z435
and indefinite permittivity D.R.S1032 . The calculations of electromagnetic
energy density in such metamaterials are of interest T.J.C205106 ; P.G.L675 ;
P.G.L075016 ; J.C163493 ; A.M634 . However, the contradictory results of
stored electromagnetic energy density in wire-SRR metamaterials have been the
subject of theoretical controversy R.R309 ; S.A.T231 ; P.M.T92 ; A.D.B165110 .
Note that in the absence of damping effects, the electromagnetic energy
density of such media can be obtained correctly by considering the
adiabatically varying electromagnetic field L.D.L . In the presence of damping
effects in wire-SRR metamaterials, on the other hand, difficulties arise in
attempting to calculate the energy associated with an electromagnetic wave
passing through such systems, and different results by two different methods
were obtained by several researchers R.R309 ; S.A.T231 ; P.M.T92 ; A.D.B165110
. The question then becomes, what is the correct result for the energy density
in wire-SRR metamaterials? In this way, one of us derived the electromagnetic
energy density formula consistent with the Landau formula in the mentioned
media P.G.L046601 .

A hyperbolic metamaterial (HMM) refers to an extremely anisotropic uniaxial
optical media and has received much attention recently C.R.S4229 ; A.P948 ;
P.S14 ; E.N214 ; L.F1 ; O.T463001 ; I.I.S . An HMM is an anisotropic medium
that has opposite signs of the two principal values of the permittivity tensor
along and perpendicular to the optic axis. Mathematically, the dispersion of
an electromagnetic wave in such a medium has a hyperbolic shape. By
definition, there are two types of HMMs. Type I with a predominantly
dielectric nature, and type II with a predominantly metallic nature. In an
appropriate frequency region, a periodic array of metallic nanowires embedded
in a host dielectric matrix is an example of an HMM of type I, and a
multilayer medium consisting of dielectric and metallic layers is an example
of an HMM of type II.

Recently, one of us P.G.L863  studied the propagation of electromagnetic
energy in the multilayer HMMs in a simple case, when metallic layers have a
permittivity of standard Drude model form, i.e., the constant part is assumed
to be one. The results show that a multilayer HMM (as an HMM of type II in the
appropriate frequency region) behaves like Drude and Lorentz medium,
respectively, when the electric field is parallel and perpendicular to the
layers. However, the problem of electromagnetic energy density in a nanowire
HMM (as an HMM of type I) have not yet been studied. This is the main
motivation of our present study. In this way, we obtain the energy density
associated with an electromagnetic wave passing through an HMM, using a simple
approach. In particular, both types of HMM are considered. For benefit of the
reader, we also show that the electromagnetic response of a composite of
metallic nanospheres embedded in a host matrix is similar to that of a Lorentz
medium (see Appendix).

##  II Basic equations

Let us consider a dispersive and absorbing HMM as an effective uniaxial
crystal. In the presence of an oscillating electric field
$\textbf{E}=\textbf{E}_{\parallel}+\textbf{E}_{\bot}$ of frequency $\omega$
and long-wavelength vibration, the motion of the effective lattice of this
effective uniaxial crystal may be represented by the Lorentz equations R.L233

$m\ddot{\textbf{r}}_{\parallel}+m\gamma\dot{\textbf{r}}_{\parallel}+m\omega_{\mathrm{0}\parallel}^{2}\textbf{r}_{\parallel}=q\textbf{E}_{\parallel}\;,$
(1)
$m\ddot{\textbf{r}}_{\bot}+m\gamma\dot{\textbf{r}}_{\bot}+m\omega_{\mathrm{0}\bot}^{2}\textbf{r}_{\bot}=q\textbf{E}_{\bot}\;,$
(2)

where $q$ is the electric charge, $m$ is the effective mass of each electric
charge, $\textbf{r}=\textbf{r}_{\parallel}+\textbf{r}_{\bot}$ is the
displacement of the oscillators, $\omega_{\mathrm{0}\bot}$
($\omega_{\mathrm{0}\parallel}$) is the resonance frequency of the charges and
$\gamma$ is the damping frequency. Note $E_{x}=E_{y}=E_{\bot}$ and
$E_{z}=E_{\parallel}$ and the $z$-axis is along the optical axis of the
effective uniaxial crystal.

Suppose that there are $N$ oscillators in the effective uniaxial crystal
volume $V$ and let us represent the effect of the high resonances by a real
constant background dielectric constants $\varepsilon_{\infty\parallel}$ and
$\varepsilon_{\infty\bot}$. Then, we have

$\textbf{P}_{\parallel}=\dfrac{Nq}{V}\textbf{r}_{\parallel}+\varepsilon_{0}\left(\varepsilon_{\infty\parallel}-1\right)\textbf{E}_{\parallel}\;,$
(3)
$\textbf{P}_{\bot}=\dfrac{Nq}{V}\textbf{r}_{\bot}+\varepsilon_{0}\left(\varepsilon_{\infty\bot}-1\right)\textbf{E}_{\bot}\;,$
(4)

and the components of the relative permittivity tensor determined by Eqs.
(1)-(4) are

$\varepsilon_{\parallel}=1+\dfrac{P_{\parallel}}{\varepsilon_{0}E_{\parallel}}=\varepsilon_{\infty\parallel}\left(1-\dfrac{F_{\parallel}\omega_{0\parallel}^{2}}{\omega(\omega+i\gamma)-\omega_{0\parallel}^{2}}\right)\;,$
(5)
$\varepsilon_{\bot}=1+\dfrac{P_{\bot}}{\varepsilon_{0}E_{\bot}}=\varepsilon_{\infty\bot}\left(1-\dfrac{F_{\bot}\omega_{0\bot}^{2}}{\omega(\omega+i\gamma)-\omega_{0\bot}^{2}}\right)\;,$
(6)

where $F_{\parallel}=\omega_{\mathrm{p}\parallel}^{2}/\omega_{0\parallel}^{2}$
($F_{\bot}=\omega_{\mathrm{p}\bot}^{2}/\omega_{0\bot}^{2}$) measures the
strength of the effective lattice resonance and
$\omega_{\mathrm{p}\parallel}=\left(Nq^{2}/\varepsilon_{0}\varepsilon_{\infty\parallel}mV\right)^{1/2}$
and
$\omega_{\mathrm{p}\bot}=\left(Nq^{2}/\varepsilon_{0}\varepsilon_{\infty\bot}mV\right)^{1/2}$.
Using the derived formula by Loudon R.L233  for the energy and damping
characteristics of the bulk electromagnetic waves in the bulk crystals, the
time-averaged total energy density associated with a wave can be written

$U=U_{\mathrm{EM}}+U_{\mathrm{K}}+U_{\mathrm{P}}\;,$ (7)

where

$U_{\mathrm{EM}}=\dfrac{1}{4}\left[\varepsilon_{0}\varepsilon_{\infty\parallel}|\textbf{E}_{\parallel}|^{2}+\varepsilon_{0}\varepsilon_{\infty\bot}|\textbf{E}_{\bot}|^{2}+\mu_{0}|\textbf{H}|^{2}\right]\;,$
(8)

is the electromagnetic energy density

$U_{\mathrm{K}}=\dfrac{1}{4}\dfrac{Nm}{V}\left(|\dot{\textbf{r}}_{\parallel}|^{2}+|\dot{\textbf{r}}_{\bot}|^{2}\right)\;,$
(9)

is the kinetic energy density of the effective uniaxial crystal, and

$U_{\mathrm{P}}=\dfrac{1}{4}\dfrac{Nm}{V}\left(\omega_{0\parallel}^{2}|\textbf{r}_{\parallel}|^{2}+\omega_{0\bot}^{2}|\textbf{r}_{\bot}|^{2}\right)\;.$
(10)

is the potential energy density of the effective uniaxial crystal. Therefore,
in general the total energy includes two parts: the energies from the electric
and magnetic fields themselves, and the energies from the medium response,
i.e., the kinetic and potential energies of the charges present under the
influence of the electromagnetic wave X.Z1853 . Also, the time-averaged power
loss density can be written as

$\displaystyle P_{\mathrm{loss}}$ $\displaystyle=$
$\displaystyle\dfrac{Nm\gamma}{2V}\left(|\dot{\textbf{r}}_{\parallel}|^{2}+|\dot{\textbf{r}}_{\bot}|^{2}\right)$
(11) $\displaystyle=$ $\displaystyle 2\gamma U_{\mathrm{K}}\;.$

Using Eqs. (1) and (2) to eliminate $\textbf{r}_{\parallel}$ and
$\textbf{r}_{\bot}$ from Eqs. (7) and (11), we obtain

$U=\dfrac{1}{4}\left[\Xi_{\parallel}|\textbf{E}_{\parallel}|^{2}+\Xi_{\bot}|\textbf{E}_{\bot}|^{2}+\mu_{0}|\textbf{H}|^{2}\right]\;,$
(12)

where

$\Xi_{\parallel}=\varepsilon_{0}\varepsilon_{\infty\parallel}\left(1+\dfrac{F_{\parallel}\omega_{0\parallel}^{2}\left(\omega^{2}+\omega_{0\parallel}^{2}\right)}{\left(\omega^{2}-\omega_{0\parallel}^{2}\right)^{2}+\gamma^{2}\omega^{2}}\right)\;,$
$\Xi_{\bot}=\varepsilon_{0}\varepsilon_{\infty\bot}\left(1+\dfrac{F_{\bot}\omega_{0\bot}^{2}\left(\omega^{2}+\omega_{0\bot}^{2}\right)}{\left(\omega^{2}-\omega_{0\bot}^{2}\right)^{2}+\gamma^{2}\omega^{2}}\right)\;,$

and

$\displaystyle P_{\mathrm{loss}}$ $\displaystyle=$
$\displaystyle\dfrac{\gamma\varepsilon_{0}}{2}\left[\Gamma_{\parallel}|\textbf{E}_{\parallel}|^{2}+\Gamma_{\bot}|\textbf{E}_{\bot}|^{2}\right]$
(13) $\displaystyle=$
$\displaystyle\dfrac{\omega\varepsilon_{0}}{2}\operatorname{Im}\left[\varepsilon_{\parallel}|\textbf{E}_{\parallel}|^{2}+\varepsilon_{\bot}|\textbf{E}_{\bot}|^{2}\right]\;,$

where

$\Gamma_{\parallel}=\varepsilon_{\infty\parallel}\dfrac{F_{\parallel}\omega_{0\parallel}^{2}\omega^{2}}{\left(\omega^{2}-\omega_{0\parallel}^{2}\right)^{2}+\gamma^{2}\omega^{2}}\;,$
$\Gamma_{\bot}=\varepsilon_{\infty\bot}\dfrac{F_{\bot}\omega_{0\bot}^{2}\omega^{2}}{\left(\omega^{2}-\omega_{0\bot}^{2}\right)^{2}+\gamma^{2}\omega^{2}}\;.$

If the losses are negligible, the time-averaged electromagnetic energy
density, for a monochromatic (single frequency) electromagnetic field, is
L.D.L

$U=\dfrac{\varepsilon_{0}}{4}\left[\dfrac{d(\omega\varepsilon_{\parallel})}{d\omega}|\textbf{E}_{\parallel}|^{2}+\dfrac{d(\omega\varepsilon_{\bot})}{d\omega}|\textbf{E}_{\bot}|^{2}\right]+\dfrac{\mu_{0}}{4}|\textbf{H}|^{2}\;.$
(14)

Setting $\gamma=0$ in Eqs. (5) and (6), and using Eq. (14), leads to

$U=\dfrac{1}{4}\varepsilon_{0}\varepsilon_{\infty\parallel}\left(1+\dfrac{F_{\parallel}\omega_{0\parallel}^{2}\left(\omega^{2}+\omega_{0\parallel}^{2}\right)}{\left(\omega^{2}-\omega_{0\parallel}^{2}\right)^{2}}\right)|\textbf{E}_{\parallel}|^{2}\\\
+\dfrac{1}{4}\varepsilon_{0}\varepsilon_{\infty\bot}\left(1+\dfrac{F_{\bot}\omega_{0\bot}^{2}\left(\omega^{2}+\omega_{0\bot}^{2}\right)}{\left(\omega^{2}-\omega_{0\bot}^{2}\right)^{2}}\right)|\textbf{E}_{\bot}|^{2}\\\
+\dfrac{1}{4}\mu_{0}|\textbf{H}|^{2}\;,$ (15)

that is equal to the result as that obtained by setting $\gamma=0$ in Eq.
(12).

##  III Nanowire hyperbolic metamaterials

Consider a periodic array of metallic nanowires with axes parallel to $z$-axis
embedded in a host dielectric matrix with the dielectric constant
$\varepsilon_{\mathrm{d}}$. The $z$-axis is along the optical axis. Let $f$ be
the filling fraction of the metallic nanowires in a unit cell satisfying
$0<f<1$. This periodic array can be used to construct an electric HMM of type
I (in an appropriate frequency region), with effective permittivity given by

$\displaystyle\varepsilon_{\mathrm{\parallel}}$ $\displaystyle=$
$\displaystyle f\varepsilon_{\mathrm{m}}+(1-f)\varepsilon_{\mathrm{d}}\;,$
(16) $\displaystyle\varepsilon_{\bot}$ $\displaystyle=$
$\displaystyle\varepsilon_{\mathrm{d}}\dfrac{(1+f)\varepsilon_{\mathrm{m}}+(1-f)\varepsilon_{\mathrm{d}}}{(1-f)\varepsilon_{\mathrm{m}}+(1+f)\varepsilon_{\mathrm{d}}}\;.$
(17)

where

$\varepsilon_{\mathrm{m}}=\varepsilon_{\infty}\left(1-\dfrac{\omega_{\mathrm{p}0}^{2}}{\omega\left(\omega+i\gamma\right)}\right)\;,$
(18)

shows the relative permittivity of a metallic nanowire with
$\varepsilon_{\infty}$ that is the permittivity in high-frequency and
$\omega_{\mathrm{p}0}$ that is the plasma frequency. Here we have relaxed the
restrictions in the previous investigation P.G.L863  by considering
$\varepsilon_{\infty}$ in Eq. (18), therefore the formulas might provide more
realistic applications.

Now, we can rewrite Eqs. (16) and (17) as Eqs. (5) and (6), respectively,
where

$\displaystyle\varepsilon_{\infty\parallel}$ $\displaystyle=$ $\displaystyle
f\varepsilon_{\infty}+(1-f)\varepsilon_{\mathrm{d}}\;,$ (19)
$\displaystyle\varepsilon_{\infty\bot}$ $\displaystyle=$
$\displaystyle\varepsilon_{\mathrm{d}}\dfrac{(1+f)\varepsilon_{\infty}+(1-f)\varepsilon_{\mathrm{d}}}{(1-f)\varepsilon_{\infty}+(1+f)\varepsilon_{\mathrm{d}}}\;,$
(20) $\displaystyle F_{\bot}$ $\displaystyle=$
$\displaystyle\dfrac{4f\varepsilon_{\mathrm{d}}}{(1-f)\left[(1+f)\varepsilon_{\infty}+(1-f)\varepsilon_{\mathrm{d}}\right]}\;,$
(21) $\displaystyle\omega_{\mathrm{p}\parallel}^{2}$ $\displaystyle=$
$\displaystyle\dfrac{f\varepsilon_{\infty}\omega_{\mathrm{p}0}^{2}}{f\varepsilon_{\infty}+(1-f)\varepsilon_{\mathrm{d}}}\;,$
(22) $\displaystyle\omega_{0\bot}^{2}$ $\displaystyle=$
$\displaystyle\dfrac{(1-f)\varepsilon_{\infty}\omega_{\mathrm{p}0}^{2}}{(1-f)\varepsilon_{\infty}+(1+f)\varepsilon_{\mathrm{d}}}\;,$
(23) $\displaystyle\omega_{0\parallel}$ $\displaystyle=$ $\displaystyle 0\;.$
(24)

These equations indicate that the responses of a nanowire HMM in the
directions perpendicular to and parallel to the optical axis are similar to
those of Lorentz and Drude media, respectively. Using the above relations, we
can obtain energy relations of HMMs of type I according to the Eqs. (7)-(15)
of the previous section. For the case $f\ll 1$, nanowires are well separated
by large distances and do not feel each other by electromagnetic interactions.
Then from Eqs. (19)-(23) we find:

$\varepsilon_{\infty\parallel}=\varepsilon_{\mathrm{d}}\;,\;\;\;\;\varepsilon_{\infty\bot}=\varepsilon_{\mathrm{d}}\;,$
$\omega_{\mathrm{p}\parallel}^{2}=\dfrac{f\varepsilon_{\infty}\omega_{\mathrm{p}0}^{2}}{\varepsilon_{\mathrm{d}}}\;,\;\;\;\;\;\omega_{0\bot}^{2}=\dfrac{\varepsilon_{\infty}\omega_{\mathrm{p}0}^{2}}{\varepsilon_{\infty}+\varepsilon_{\mathrm{d}}}\;,\;\;\;\;\;F_{\bot}=\dfrac{4f\varepsilon_{\mathrm{d}}}{\varepsilon_{\infty}+\varepsilon_{\mathrm{d}}}\;,$
(25)

where
$\varepsilon_{\infty}\omega_{\mathrm{p}0}/\sqrt{\varepsilon_{\infty}+\varepsilon_{\mathrm{d}}}$
is the frequency of dipolar resonance of a single metallic nanowire surrounded
by a dielectric medium A.M064502 .

Note that for an electric HMM of type I we should have ${\rm
Re}\left[\varepsilon_{\mathrm{\parallel}}\right]<0$ and ${\rm
Re}\left[\varepsilon_{\mathrm{\bot}}\right]>0$ and the imaginary parts of them
must be small enough to be negligible. Therefore, using Eqs. (5) and (6) we
find ${\rm Re}\left[\varepsilon_{\mathrm{\parallel}}\right]<0$, if
$\omega<\sqrt{\omega_{\mathrm{p}\parallel}^{2}-\gamma^{2}}$, and ${\rm
Re}\left[\varepsilon_{\mathrm{\bot}}\right]>0$, if $\omega<\omega_{0\bot}$,
that means

$\omega<\mathrm{min}\left(\sqrt{\omega_{\mathrm{p}\parallel}^{2}-\gamma^{2}},\omega_{0\bot}\right)\;.$
(26)

##  IV Multilayer hyperbolic metamaterials

Consider a multilayer structure consisting of isotropic metal/dielectric
layers. The $z$-axis is directed along the optical axis. Let
$\varepsilon_{\mathrm{m}}$ ($\varepsilon_{\mathrm{d}}$) be the relative
permittivity of the metal (dielectric) layer, and let $f$ be the filling ratio
of the metal layer satisfying $0<f<1$. Also, Eq. (18) shows the relative
permittivity of a metallic layer. This multilayer structure can be used to
construct an electric HMM of type II (in an appropriate frequency region),
with effective permittivity given by

$\displaystyle\varepsilon_{\mathrm{\parallel}}$ $\displaystyle=$
$\displaystyle\dfrac{\varepsilon_{\mathrm{m}}\varepsilon_{\mathrm{d}}}{(1-f)\varepsilon_{\mathrm{m}}+f\varepsilon_{\mathrm{d}}}\;,$
(27) $\displaystyle\varepsilon_{\mathrm{\bot}}$ $\displaystyle=$
$\displaystyle f\varepsilon_{\mathrm{m}}+(1-f)\varepsilon_{\mathrm{d}}\;.$
(28)

We can rewrite Eqs. (27) and (28) as Eqs. (5) and (6), respectively, where
P.G.L863

$\displaystyle\varepsilon_{\infty\bot}$ $\displaystyle=$ $\displaystyle
f\varepsilon_{\infty}+(1-f)\varepsilon_{\mathrm{d}}\;,$ (29)
$\displaystyle\varepsilon_{\infty\parallel}$ $\displaystyle=$
$\displaystyle\dfrac{\varepsilon_{\infty}\varepsilon_{\mathrm{d}}}{(1-f)\varepsilon_{\infty}+f\varepsilon_{\mathrm{d}}}\;,$
(30) $\displaystyle F_{\parallel}$ $\displaystyle=$
$\displaystyle\dfrac{f}{1-f}\dfrac{\varepsilon_{\mathrm{d}}}{\varepsilon_{\infty}}\;,$
(31) $\displaystyle\omega_{\mathrm{p}\bot}^{2}$ $\displaystyle=$
$\displaystyle\dfrac{f\varepsilon_{\infty}\omega_{\mathrm{p}0}^{2}}{f\varepsilon_{\infty}+(1-f)\varepsilon_{\mathrm{d}}}\;,$
(32) $\displaystyle\omega_{0\parallel}^{2}$ $\displaystyle=$
$\displaystyle\dfrac{(1-f)\varepsilon_{\infty}\omega_{\mathrm{p}0}^{2}}{(1-f)\varepsilon_{\infty}+f\varepsilon_{\mathrm{d}}}\;,$
(33) $\displaystyle\omega_{0\bot}$ $\displaystyle=$ $\displaystyle 0\;.$ (34)

These equations indicate that the responses of a multilayer HMM in the
directions perpendicular to and parallel to the optical axis are similar to
those of Drude and Lorentz media, respectively. Using the above relations, we
can obtain energy relations of HMMs of type II according to the Eqs. (7)-(15)
of Sec. II. Again, if we consider $f\ll 1$ then from Eqs. (29)-(33) we find:

$\varepsilon_{\infty\bot}=\varepsilon_{\mathrm{d}}\;,\;\;\;\;\varepsilon_{\infty\parallel}=\varepsilon_{\mathrm{d}}\;,$
$\omega_{\mathrm{p}\bot}^{2}=\dfrac{f\varepsilon_{\infty}\omega_{\mathrm{p}0}^{2}}{\varepsilon_{\mathrm{d}}}\;,\;\;\;\;\;\omega_{0\parallel}^{2}=\omega_{\mathrm{p}0}^{2}\;,\;\;\;\;\;F_{\parallel}=f\dfrac{\varepsilon_{\mathrm{d}}}{\varepsilon_{\infty}}\;,$
(35)

where $\omega_{\mathrm{p}0}$ is the resonance frequency of a single metallic
layer in a dielectric medium. This indicates that the resonance effect also
happens in the metallic layers of the system, however as a result we note that
this resonance frequency is free from the effect of permittivity of the
dielectric medium.

Figure 1: Variation of $\mathcal{D}_{\mathrm{energy}}$ (blue curves), and
$\mathcal{D}_{\mathrm{loss}}$ (red curves) with respect to the dimensionless
frequency $\omega/\omega_{\mathrm{p}0}$ for a nanowire HMM as a type I HHM.
Here, the filling fraction is $f=0.2$, the relative permittivity of the host
dielectric matrix is $\varepsilon_{\mathrm{d}}=2.5$, $\varepsilon_{\infty}=1$,
and $\gamma=0.01\omega_{\mathrm{p}0}$. Therefore, we have
$\varepsilon_{\infty\parallel}=2.2$, $\varepsilon_{\infty\bot}=2.11$,
$\omega_{\mathrm{p}\parallel}=0.3\omega_{\mathrm{p}0}$,
$\omega_{0\bot}=0.459\omega_{\mathrm{p}0}$, and $F_{\bot}=0.78$. Also, the
vertical dashed line shows $\omega/\omega_{\mathrm{p}0}=0.301$. Figure 2:
Variation of $\mathcal{D}_{\mathrm{energy}}$ (blue curves), and
$\mathcal{D}_{\mathrm{loss}}$ (red curves) with respect to the dimensionless
frequency $\omega/\omega_{\mathrm{p}0}$ for a multilayer HMM as a type II HHM.
Here, the filling fraction is $f=0.2$, the relative permittivity of the
dielectric layers is $\varepsilon_{\mathrm{d}}=2.5$, $\varepsilon_{\infty}=1$,
and $\gamma=0.01\omega_{\mathrm{p}0}$. Therefore, we have
$\varepsilon_{\infty\bot}=2.2$, $\varepsilon_{\infty\parallel}=1.92$,
$\omega_{\mathrm{p}\bot}=0.3\omega_{\mathrm{p}0}$,
$\omega_{0\parallel}=0.78\omega_{\mathrm{p}0}$, and $F_{\parallel}=0.625$.
Also, the vertical dashed line shows $\omega/\omega_{\mathrm{p}0}=0.301$.

Note that for an electric HMM of type II we should have ${\rm
Re}\left[\varepsilon_{\mathrm{\parallel}}\right]>0$ and ${\rm
Re}\left[\varepsilon_{\mathrm{\bot}}\right]<0$. Therefore, using Eqs. (5) and
(6) we find ${\rm Re}\left[\varepsilon_{\parallel}\right]>0$, if
$\omega<\omega_{0\parallel}$, and ${\rm Re}\left[\varepsilon_{\bot}\right]<0$,
if $\omega<\sqrt{\omega_{\mathrm{p}\bot}^{2}-\gamma^{2}}$, that means

$\omega<\mathrm{min}\left(\sqrt{\omega_{\mathrm{p}\bot}^{2}-\gamma^{2}},\omega_{0\parallel}\right)\;.$
(36)

##  V Numerical results

In this section, numerical results for the distribution of electromagnetic
energy density in an HMM are presented. We assume that the filling fraction is
$f=0.2$, the relative permittivity of the dielectric medium is
$\varepsilon_{\mathrm{d}}=2.5$, $\varepsilon_{\infty}=1$, and
$\gamma=0.01\omega_{\mathrm{p}0}$. Therefore, we use the following values for
the parameters appearing in Eqs. (12) and (13):
$\varepsilon_{\infty\parallel}=2.2$, $\varepsilon_{\infty\bot}=2.11$,
$\omega_{\mathrm{p}\parallel}=0.3\omega_{\mathrm{p}0}$,
$\omega_{0\bot}=0.459\omega_{\mathrm{p}0}$, and $F_{\bot}=0.78$ for a periodic
array of metallic nanowires as a type I HMM, and
$\varepsilon_{\infty\bot}=2.2$, $\varepsilon_{\infty\parallel}=1.92$
$\omega_{\mathrm{p}\bot}=0.3\omega_{\mathrm{p}0}$,
$\omega_{0\parallel}=0.78\omega_{\mathrm{p}0}$, and $F_{\parallel}=0.625$ for
a multilayer HMM as a type II HMM.

To study the corresponding dispersion/anisotropy and absorption effects in the
energy density and power loss, we use the directionality factors P.G.L863 , as
$\mathcal{D}_{\mathrm{energy}}=\Xi_{\bot}/\Xi_{\parallel}$, where
$\Xi_{\parallel}$ and $\Xi_{\bot}$ are the coefficients
$|\textbf{E}_{\parallel}|^{2}$ and $|\textbf{E}_{\bot}|^{2}$ in Eq. (12),
respectively, and
$\mathcal{D}_{\mathrm{loss}}=\Gamma_{\bot}/\Gamma_{\parallel}$, where
$\Gamma_{\parallel}$ and $\Gamma_{\bot}$ are the coefficients
$|\textbf{E}_{\parallel}|^{2}$ and $|\textbf{E}_{\bot}|^{2}$ in Eq. (13),
respectively.

Figs. 1 and 2 show the variation of $\mathcal{D}_{\mathrm{energy}}$ (blue
curves), and $\mathcal{D}_{\mathrm{loss}}$ (red curves) with respect to the
dimensionless frequency $\omega/\omega_{\mathrm{p}0}$ for HHMs of types I and
II, respectively. From Fig. 1 one can see the absorption property of type I
HMM becomes isotropic with respect to the electric field at

$\dfrac{\omega}{\omega_{\mathrm{p}0}}=\sqrt{\dfrac{f^{2}}{\left[f+(1-f)\varepsilon_{\mathrm{d}}\right]^{2}}-\gamma^{2}}=0.301\;,$

while the energy storage ability of the medium becomes isotropic with respect
to the electric field for $\omega/\omega_{\mathrm{p}0}<0.301$. Also,
$\mathcal{D}_{\mathrm{loss}}$ and $\mathcal{D}_{\mathrm{energy}}$ decrease
with decreasing the value of $\omega/\omega_{\mathrm{p}0}$. Furthermore, from
Fig. 2 it is clear that the energy storage ability of the type II HMM becomes
almost isotropic with respect to the electric field at
$\omega/\omega_{\mathrm{p}0}=0.301$. Below this frequency, the type II HMM can
store more energy if the electric field direction is parallel to the layers.

##  VI Conclusion

Although it seems that the structure of an HMM is much simpler than the wire-
SRR and chiral metamaterials, however, the effective permittivities of HMMs
are more complex than those of wire-SRR and chiral metamaterials and
investigation of electromagnetic energy density in such media may be
difficult. To remedy this difficulty, we have obtained familiar forms for the
effective permittivities of nanowire HMMs and multilayer HMMs, i.e., similar
to the Lorentz and Drude media. These new forms for the effective
permittivities of the HMMs simply show that HMMs have different dynamical
properties in the direction parallel and perpendicular to their optical axis.
Using these new forms of the effective permittivities of HMMs, we have derived
the energy density associated with an electromagnetic wave passing through a
HMM. In this way, HMM of types I and II have been treated, separately.

## Appendix: Lorentz model of a composite of metallic nanospheres

Consider a composite of Drude metallic nanospheres with relative dielectric
constant
$\varepsilon_{\mathrm{m}}=\varepsilon_{\infty}\left(1-\omega_{\mathrm{p}0}^{2}/\omega\left(\omega+i\gamma\right)\right)$
embedded in a host matrix with relative dielectric constant
$\varepsilon_{\mathrm{d}}$. Let $f$ be the volume fraction of the embedded
nanospheres satisfying $0<f<1$. For the relative effective permittivity of the
present composite we have U.K ; A.M

$\dfrac{\varepsilon_{\mathrm{eff}}-\varepsilon_{\mathrm{d}}}{\varepsilon_{\mathrm{eff}}+2\varepsilon_{\mathrm{d}}}=f\dfrac{\varepsilon_{\mathrm{m}}-\varepsilon_{\mathrm{d}}}{\varepsilon_{\mathrm{m}}+2\varepsilon_{\mathrm{d}}}\;,$
(A-1)

or

$\varepsilon_{\mathrm{eff}}=\varepsilon_{\mathrm{d}}\dfrac{\left(\varepsilon_{\mathrm{m}}+2\varepsilon_{\mathrm{d}}\right)+2f\left(\varepsilon_{\mathrm{m}}-\varepsilon_{\mathrm{d}}\right)}{\left(\varepsilon_{\mathrm{m}}+2\varepsilon_{\mathrm{d}}\right)-f\left(\varepsilon_{\mathrm{m}}-\varepsilon_{\mathrm{d}}\right)}\;.$
(A-2)

We can rewrite Eq. (A-2) as

$\varepsilon_{\mathrm{eff}}=\varepsilon_{\mathrm{b}}\left(1-\dfrac{F\omega_{0}^{2}}{\omega(\omega+i\gamma)-\omega_{0}^{2}}\right)\;,$
(A-3)

where

$\displaystyle\varepsilon_{\mathrm{b}}$ $\displaystyle=$
$\displaystyle\varepsilon_{\mathrm{d}}\dfrac{\left(\varepsilon_{\infty}+2\varepsilon_{\mathrm{d}}\right)+2f\left(\varepsilon_{\infty}-\varepsilon_{\mathrm{d}}\right)}{\left(\varepsilon_{\infty}+2\varepsilon_{\mathrm{d}}\right)-f\left(\varepsilon_{\infty}-\varepsilon_{\mathrm{d}}\right)}\;,$
(A-4) $\displaystyle F$ $\displaystyle=$
$\displaystyle\dfrac{9f\varepsilon_{\mathrm{d}}}{(1-f)\left[\left(\varepsilon_{\infty}+2\varepsilon_{\mathrm{d}}\right)+2f\left(\varepsilon_{\infty}-\varepsilon_{\mathrm{d}}\right)\right]}\;,$
(A-5) $\displaystyle\omega_{0}^{2}$ $\displaystyle=$
$\displaystyle\dfrac{(1-f)\varepsilon_{\infty}\omega_{\mathrm{p}0}^{2}}{\left(\varepsilon_{\infty}+2\varepsilon_{\mathrm{d}}\right)-f\left(\varepsilon_{\infty}-\varepsilon_{\mathrm{d}}\right)}\;.$
(A-6)

For the case $f\ll 1$ we find
$\omega_{0}=\varepsilon_{\infty}\omega_{\mathrm{p}0}/\sqrt{\varepsilon_{\infty}+2\varepsilon_{\mathrm{d}}}$
that is the frequency of dipolar resonance of a single metallic nanosphere
surrounded by a dielectric medium A.M849 .

## References

  * (1) D. R. Smith, W. J. Padilla, D. C. Vier, S. C. Nemat-Nasser, S. Schultz, Phys. Rev. Lett. 84, 4184 (2000)
  * (2) J. B. Pendry, Phys. Rev. Lett. 85, 3966 (2000)
  * (3) X. Zhang, Z. Liu, Nat. Mater. 7, 435 (2008)
  * (4) D. R. Smith, P. Kolinko, D. Schurig, J. Opt. Soc. Am. B 21, 1032 (2004)
  * (5) T. J. Cui, J. A. Kong, Phys. Rev. B 70, 205106 (2004)
  * (6) P. G. Luan, Y. T. Wang, S. Zhang, X. Zhang, Opt. Lett. 36, 675 (2011)
  * (7) P. G. Luan, J. Phys. Commun. 2, 075016 (2018)
  * (8) J. Chen, Y. Xuanyuan, Y. Dai, Optik, 207, 163493 (2020)
  * (9) A. Moradi, J. Mod. Opt. 68, 634 (2021)
  * (10) R. Ruppin, Phys. Lett. A 299, 309 (2002)
  * (11) S. A. Tretyakov, Phys. Lett. A 343, 231 (2005)
  * (12) P. M. T. Ikonen, S. A. Tretyakov, IEEE Trans. Microwave Theory Tech. 55, 92 (2007)
  * (13) A. D. Boardman, K. Marinov, Phys. Rev. B 73, 165110 (2006)
  * (14) L. D. Landau, E. M. Lifshitz, Electrodynamics of Continuous Media, 2nd ed. (Pergamon Press, New York, 1984)
  * (15) P. G. Luan, Phys. Rev. E 80, 046601 (2009)
  * (16) C. R. Simovski, P. A. Belov, A. V. Atrashchenko, Y. S. Kivshar, Adv. Mater. 24, 4229 (2012)
  * (17) A. Poddubny, I. Iorsh, P. Belov, Y. Kivshar, Nature Photon 7, 948 (2013)
  * (18) P. Shekhar, J. Atkinson, Z. Jacob, Nano Convergence 1, 14 (2014)
  * (19) E. Narimanov, A. Kildishev, Nature Photon 9, 214 (2015)
  * (20) L. Ferrari, C. Wu, D. Lepage, X. Zhang, Z. Liu, Prog. Quantum Electron 40, 1 (2015)
  * (21) O. Takayama, A. A. Bogdanov, A. V. Lavrinenko, J. Phys.: Condens. Matter 29, 463001 (2017)
  * (22) I. I. Smolyaninov, Hyperbolic Metamaterials (Morgan and Claypool/Institute of Physics, London, 2018)
  * (23) P.-G. Luan, Crystals 10, 863 (2020)
  * (24) R. Loudon, J. Phys. A: Gen. Phys. 3, 233 (1970)
  * (25) X. Zheng, P. Palffy-Muhoray, Phys. Lett. A 379, 1853 (2015)
  * (26) A. Moradi, Phys. Plasmas 27, 064502 (2020)
  * (27) U. Kreibig, M. Vollmer, Optical Properties of Metal Clusters, (Springer, Berlin, 1995)
  * (28) A. Moradi, Canonical Problems in the Theory of Plasmonics: From 3D to 2D Systems, (Springer, Switzerland, 2020)
  * (29) A. Moradi, Z. Naturforsch. 71, 849 (2016)