VO$_2$ as a natural optical metamaterial

VO$_2$ is a unique phase change material with strongly anisotropic electronic properties. Recently, samples have been prepared that present a co-existence of phases and thus form metal-insulator junctions of the same chemical compound. Using first principles calculations, the optical properties of metallic and semiconducting VO$_2$ are here discussed to design self-contained natural optical metamaterials, avoiding coupling with other dielectric media. The analysis of the optical properties complements the experiments in the description of the vast change in reflectance and metallicity for both disordered and planar compounds. The present results also predict the possibility to realize ordered VO$_2$ junctions operating as efficient hyperbolic metamaterials in the THz-visible range, by simply adjusting the ratio between metallic and insulating VO$_2$ content. The possibility to excite propagating {\em volume plasmom polariton} across the metamaterial is finally discussed.

VO2 is a unique phase change material with strongly anisotropic electronic properties. Recently, samples have been prepared that present a co-existence of phases and thus form metal-insulator junctions of the same chemical compound. Using first principles calculations, the optical properties of metallic and semiconducting VO2 are here discussed to design self-contained natural optical metamaterials, avoiding coupling with other dielectric media. The analysis of the optical properties complements the experiments in the description of the vast change in reflectance and metallicity for both disordered and planar compounds. The present results also predict the possibility to realize ordered VO2 junctions operating as efficient hyperbolic metamaterials in the THz-visible range, by simply adjusting the ratio between metallic and insulating VO2 content. The possibility to excite propagating volume plasmom polariton across the metamaterial is finally discussed.

I. INTRODUCTION
Optical metamaterials (OMMs) 1 are artificially nanostructured compounds that exhibit an unconventional optical response to light that is unavailable in natural materials. Although the field of metamaterials started with the search of the so-called double negative materials (or Veselago media) 2 , the idea of realizing artificial materials that can be engineered to support unique electromagnetic modes in the subwavelength regime opened up new conceptual frontiers for OMMs, which rapidly surpassed the concept of negative refraction. This has fostered the research of a new class of tunable and active metamaterials, known as hyperbolic metamaterials (HMMs) 3 .
HMMS are highly anisotropic uniaxial materials that derive their name from the topology of the isofrequency surface ω(k)=constant. In uniaxial materials, the dielec-tricˆ and the permeabilityμ tensors can be written in a diagonal form, in terms of the four fundamental parameters ( , ⊥ , µ , µ ⊥ ), where the subscripts and ⊥ indicate components parallel and perpendicular to the anisotropy (z) axis. If one out of the four fundamental parameters is negative, the isofrequency surface opens into a hyperboloid. For non-magnetic materials, the choice x,y = ⊥ > 0, z = < 0 corresponds to a twofold hyperboloid and the medium is called a type-I metamaterial; the choice x,y = ⊥ < 0, z = > 0 describes a one-fold hyperboloid and the medium is called a type-II material. These unusual characteristics open the way to a wide range of unprecedented applications such as negative refraction 4 , sub-diffraction imaging 5 , subwavelength waveguides 6 , and spontaneous emission engineering 7 .
Hyperbolic electric dispersion is typically realized by using artificial nanostructures, such as multilayers of alternating sub-wavelength layers of metals and dielectrics or assemblies of metallic nanowires embedded in a dielectric host. The choice of the structural parameters (e.g., thickness and number of layers) as well as of the constituting materials affects the optical and plasmonic properties of the HMM as well as the frequency range of their working application. On one hand, manufacturing of artificial HMMs is still a challenge as it requires ordered growth in the sub-wavelength scale of different materials with potentially different crystalline characteristics and/or growth conditions. On the other hand, the discovery of natural HMMs 8,9 -i.e. uniaxial single-phase materials with opposite sign of permittivity -would offer a viable alternative to artificial metal-dielectric composites. However, the hyperbolic requirement ⊥ < 0 is a rare condition and only few systems (such as graphitelike materials, cuprate and ruthenate perovskites) have been identified as natural HMMs up to now 10 .
In this regard, phase change materials (PCMs) may constitute a valuable compromise for the realization of HMMs, playng with different phases, rather than different materials, thus reducing defect and strain content. Vanadium dioxide (VO 2 ) in particular has shown promise due to its metal-to-insulator transition (MIT) 11 at the critical temperature of 340K 12 . The MIT in VO 2 occurs when the atoms in the metallic structure dimerize and break the crystal symmetry, resulting in a shift from a tetragonal rutile (R) phase to a low temperature monoclinic (M 1 ) phase (Fig. 1). The optical properties differ widely between phases, making it possible to realize switchable devices based on altering the absorption, conductivity, or reflectance of the material by precise control of the system environment. In fact, schemes are already in place to implement VO 2 in switching plasmonic nanoantenna, 13 ultrafast light emission modulators 14 , near-field thermal transfer devices 15 , and novel optical metamaterials 16,17 , in connection with other dielectric oxides (e.g. sapphire 16 , titania 18 , silica 19 ) or conducting media (e.g. Au 20 , Al:ZnO 21 ).
Additionally, recent experiments have shown that it is possible to create VO 2 samples containing a co-existent mixture of both phases, i.e. to create a compound that combines the original metallic and dielectric behavior of the separate phases into a brand new metamaterial, with intrinsic optical properties different from the ones of the original constituents. This characteristic renders VO 2 a unique PCM. These VO 2 mixtures have often been induced through chemical doping 22 , temperature gradients 13,23 , strain 24 , charge injection 25 , or ion beam irradiation 16 . While several experiments have been able to accurately model the optical constants of inhomogeneous 26,27 or two-dimensional 16 VO 2 metamaterial at a target wavelength, a comprehensive study of the effective optical properties of mixed-phase VO 2 across a broad range of energy and composition is still lacking.
On the theoretical side, the majority of first principles investigations remained focused on describing the electronic band structure and characterizing the mechanism behind MIT 28-30 , while modelling of their plasmonic properties uses empiric parameters extrapolated from experiments 19,31 . In particular, the use of empirical data suffers from lack of transferability, as they are strictly dependent to the specific characteristics (i.e. structural configuration, temperature, growth environment, defects, dopants, etc) of the very sample they are extracted from. Thus a fully ab initio description of the optical properties of both VO 2 phases at the same level of accuracy would be a dramatic improvement in the understanding of these PCM materials. First-principles simulations of VO 2 optical response is a dramatic challenge because of the high electron-electron correlation deriving from partially occupied and localized V (3d) orbitals. This is particularly critical for the description of the monoclinc phase where the presence of the energy bandgap and of excitonic ef-fects have to be carefully taken into account. Different approaches have been proposed so far to treat electron correlation and optical properties of VO 2 , including dynamical mean field theory (DMFT) 28 and many-boby GW Bethe-Salpheter approaches [32][33][34] . The inclusion of quasi-particle and local-field effects have been demonstrated to correctly reproduce the low-energy spectrum and band-gap onset of the semiconducting phase. On the other hand, such advanced many-body techniques do not give any remarkable improvement to the description of the optical properties of metallic phase, with respect to computationally simpler single-particle approaches.
Here, we provide an ab initio investigation to accurately describe the optical properties of both metallic and semi-conducting VO 2 , and we then present a description of a combined metamaterial with a varying degree of phase proportions and structural arrangements. First principles calculations represent a predictive parameterfree instrument for the microscopic understanding of complex systems beyond a specific experimental condition, like applied stress, doping, defects, etc. Thus, when implementing an effective medium approximation (EMA) on the basis of the ab initio results, we are able to reproduce the vast change in reflectance and metallicity in both disordered and planar compounds, observed in the experiments 16,27,35 . Finally, we investigate the hyperbolic behavior of vertically stacked VO 2 herostructures and we characterize the angular behavior of so called volume plasmon polaritons 36 , which may propagate across the metamaterial.

II. METHOD
We employ first principles DFT approaches as implemented in the Quantum-Espresso 37 software package to calculate bulk material properties for each phase of VO 2 . The exchange-correlational functional is approached through a local density approximation (LDA), and norm-conserving pseudopotentials are used to describe ion cores. Following previous theoretical investigations [28][29][30]38 , we use experimentally obtained lattice parameters 12 for both systems and optimize internal atomic positions such that atomic forces are less than 0.01 eV/Å. Single particle wavefunctions are expanded in a planewave basis set up to an energy cut-off of 140 Ry, and initial calculations use uniform Monkorst-Pack grids of (8 × 8 × 8) and (6 × 6 × 6) for the respective R and M 1 unit cells. We additionally employ a recent pseudohybrid Hubbard implementation of DFT+U, namely ACBN0 39 that profitably corrects the DFT energy bandgap 40 as well as the dielectric and vibrational properties of metal-oxides 41 . The U values for Vanadium resulting from ACBN0 cycle and used for the calculations are U(V d ) =0.95 eV and U(V d )= 2.0 eV, for M 1 and R phases, respectively. In the monoclinic phase, two sets of inequivalent oxygen are present, whose U(O p ) values are 6.59 eV and 6.97 eV, while U(O p )= 6.93 eV is the value for oxygen in the tetragonal phase.
In the spirit of investigating the optical properties of mixed compounds, where the presence of screening metallic components make correlation effects not crucial for the overall metamaterial, the optical properties are simulated at the single-particle level starting from the DFT+U description of the electronic ground state. This assures a reasonable compromise between numerical accuracy and computational cost for both M 1 and R phase and opens to the future possibility of treating larger simulations cells able to include dopants, defects or local disordered elements. In particular, the full complex dielectric functionˆ = r + i i and loss function L = −Im[ˆ −1 (ω)] are determined using the inbuilt epsilon.x code, contained within the Quantum-Espresso package, which implements an independent-particle formulation of the frequency-dependent Drude-Lorentz model for the dielectric functionˆ (ω), where both intraband (Drude-like) and interband (Lorentz-like) contributions are explicitly considered 42 . This approach accounts only for the electronic contribution to the dielectric function: at low energies the interaction with optical phonon modes may affect the optical response of the system and, in a few cases, form phonon-polariton excitations not considered here.
Once the full complex dielectric functions for both material phases are obtained, we employ an effective medium approximation (EMA) 43 to simulate the optical constants˜ =˜ r + i˜ i of the joint-phase metamaterial. Finally, the reflectance function R of the composite is easily obtained as: are the real (ñ) and imaginary (k) part of the effective refractive index.

A. Single phase characterization
As a first step, we calculate the electronic structure of VO 2 for both the rutile and monoclinic phases, whose atomic structures are shown in Fig. 1 and bandstructure plots are displayed in Fig. 2. The R phase [panel (a)] is metallic and characterized by a Fermi level (E F ) that overlaps a manyfold of bands with a predominant V 3d character. In the octahedral coordination symmetry of rutile crystal, the V 3d states are split into σ-bonded (e σ g ) orbitals, more strongly hybridized with O, and πbonded (t 2g ) orbitals. The t 2g bands can be separated into one a 1g band crossing the Fermi level, which derives from the overlapping of V 3d orbitals along the rutile c-axis, and the remaining e π g subbands located just above E F . The lower-energy bands between -2 and -7 eV have a prominent O 2p character, only partially mixed with the V states, and are separated from t 2g bands by ∼1 eV. This matches the experimental values obtained by photoemission spectra 44,45 and previous theoretical simulations 28,30,38 .
In the monoclinic phase, the crystalline distortion from the tetragonal phase induces an alternate V-V dimerization along the c-axis [ Fig. 1(b)], in agreement with the experimental findings 11,27 . The structural symmetrybreaking reflects on the bandstructure [ Fig. 2(b)], which exhibits a negative energy shift of the lowest hybridized 3d states to just below E F , resulting in a semiconductor with a narrow indirect bandgap of 0.45 eV between Z and Γ, in agreement with previous experimental 27, 44,46 and theoretical values 28,47 . More specifically: due to the V-V dimerization the a 1g band splits into bondingantibonding subbands. This causes the lowering of the bonding-a 1g band and the upshift of the e π g subbands, opening the gap at the Fermi energy. A multiplet of lower energy O 2p -derived bands appears at ∼-2 eV, similar to the R case. Notably, for both phases the inclusion of the Hubbard-like correction pushes the upper V 3d and the O 2p bands towards higher and lower energy, respectively, leaving the t 2g bands close to Fermi energy almost unvaried. This behavior well agrees with the Mott-Peierls character of the M 1 phase 27,28,48 . VO 2 is highly anisotropic and this reflects on different optical response in the directions parallel or perpendicular to the c-axis of the structure (z direction in Fig. 1). Both parallel and perpendicular components are shown in Fig. 3. The rutile system is characteristic of a Drudelike picture: the divergence at zero energy and the rapid decay of the imaginary part of the dielectric function are representative of the energy damping associated with the dc electrical resistivity of metals. Accordingly, r starts strongly negative at low energy and becomes positive at the crossover energy E p = 2.4 eV. This, along with the condition i ≈ 0, results in a peak in the loss function (green curves), which can be directly compared with experiments. The simulated crossover energy well reproduces the experimental value of 2.75 eV detected in EELS measurements 49 . The peak at E p = 2.4 eV is interpreted as a screened plasmon excitation, typical of transition metals, including Ag and Au 50 . The origin of this plasmonic peak is associated to the coexistence of intra-and interband transitions [see, e.g., green and red arrows in Fig. 2(a)]. Free electron oscillations stem from the intraband transitions of a 1g bands that crosses the Fermi level, giving rise to the Drude-tail of r for E→ 0. At higher excitation energies, the activation of interband transitions between, e.g., occupied O(2p) and empty V (t 2g ) states gives a positive contribution to the (negative) Drude component of r that becomes positive at the crossover energy E p . This plasmon is usually labeled as screened as it involves the collective oscillation only of a fraction of the total valence electrons that can be considered as free, the rest being effectively screened by interband transitions. The broad and structured band in i spectrum for E> 2.5 eV is associated with the optical absorption of the material in the visible range and is due to interband transitions between occupied O 2p and unoccupied V 3d states. The crossover energy E p thus identifies two well distinguished optical characters for the material: reflective for E < E p and absorptive for E > E p .
The dielectric function of monoclinic phase has the behavior typical of semiconductors [ Fig. 3(b)]: the real part of the dielectric function has finite value at E = 0 eV, which corresponds to the (parallel and perpendicular) value of the static dielectric constant. The alternated V-V dimerization along the c-axis enhances the differences between the parallel and perpendicular components of the dielectric function. While ⊥ remains always positive over the entire considered energy range, changes sign in the visible range, exhibiting a typical Lorentz behavior. This generates the appearance of a crossover frequency at E=2.66 eV, which corresponds to a peak in the corresponding loss spectrum. The imaginary parts are dominated by a double peak between 1-2 eV. The lowest energy component corresponds to the onset of the insulating gap. The second major peak at 2.1 eV is the result of larger energy transitions between the O(2p) states and the lowest conduction band at Γ point.

B. Joint-phase mixtures
Although the optical anisotropy is a general condition in uniaxial systems, this is often not sufficient to qualify the intrinsic material as a natural HMM. Following the pioneering work by Qazilbash and coworkers 27,35,49 , we considered mixed phase films whose relative semiconductor/metal ratio can be reversibly modified. In particular, starting from room temperature semiconducting system they observed the appearance of metallic puddles at the critical temperature. Puddles enlarge and coalesce as the temperature is further increased to form a complete metallic phase. Provided the incident wavelength is large compared to the size of the metallic inhomogeneities, the mixed system can be described by an effective dielectric function within the EMA.
For a homogeneous disordered systems [ Fig. 4(a)] we assumed a Bruggeman mixing model (BMM) 43 , where the effective dielectric function˜ is given by solving the where m and d are the complex dielectric functions of the metal rutile and dielectric monoclinic phase, respectively, calculated from first principles 51 . In this case, when calculating˜ , the different vector components of the dielectric tensor for each individual R and M 1 phase are first averaged over the three spatial directions so as to simulate the random crystal orientations found in a realistic experimental environment. In Eq. (3) the volume fractions of the metallic and insulating phases are given by f and (1 − f ), respectively; q is a depolarization factor that depends on the geometry of the phases. If we consider an insulating film with evenly spaced, spherical metallic inclusions, the depolarization factor, q, is given by 1/3. This is a valid assumption as long as the diameter of each metallic sphere does not exceed the thickness of the film. The simulated 2D plot of reflectance R(E, f ) of the mixed phase metamaterial, as a function of the incoming radiation energy (E) and the filling fraction (f ) is shown in Fig. 5(a). For f = 0 and f = 1, R recovers the semiconducting and metallic behavior of the M 1 and R phases, respectively. The most striking feature of the reflectance spectrum is the behavior in an energy window of 0.2 -1.0 eV. If we focus on a single energy value in this range, we witness the reflectance of the sample transition from perfectly absorbing to reflecting as the percentage of metal is increased. The modulation of R is in agreement with experiments that exhibit orders of magnitude shifts in reflectance of a VO 2 sample in the near infra-red energy regime as the material undergoes the MIT 16,17,26,27 . However, in the comparison with experiment we have to notice that our results are relative to pure VO 2 films. All experiments to-date necessitate an underlying substrate that may influence the overall optical properties of the whole VO 2 /substrate system. Common substrate, such as TiO 2 or sapphire, are themselves strongly absorbing at singular wavelengths in THz region. This may explain why experiments have thus far seen near-perfect absorp- tion at only a specific energy value and not in a larger energy range.
For E < 1.0 eV, while there is a point of minimal reflectance in the medium-near IR region for each filling factor, the energy value at which the film demonstrates lowest reflectance is blue-shifted as the metallic component is increased. For f > 0.4 the system acts as an effective metallic compound. In the energy range 1.0-2.5 eV the mixture is only partially reflecting for every semiconductor/metal composition with R being 30-50% on average. Additionally, there is a point of even lower reflectance for all phase proportions in the visible spectrum, at approximately 2.25 eV. At higher energies, i.e. beyond the crossover frequency of both M 1 and R phases, the metamaterial is no more reflecting and the absorption processes become preponderant. Occurrence of metallicity at f = 0.4 is confirmed by the plot of the real and imaginary part of the effective dielectric function, shown in Fig. 5(b) for the same f value.
The spontaneous arrangement of the metallic puddles in mixed films does not allow to define preferential directions for the effective optical response. However, throughout the spatial control of defects or strain it is possible to fabricate joint-phase metamaterials with periodic sequence of metallic and semiconducting sections. In particular, following the experimental work by Rensberg et al (Ref. [ 16]), we considered a ridged heterostructure (RHS) of alternating stripes of metallic and semiconducting VO 2 material, as shown in Fig. 4b. In this case, we can easily distinguish two components of the dielectric The structural spatial anisotropy of ridged heterostructure clearly reflects on the optical response, as shown in Fig. 6. At low energies (E< 0.5 eV), when polarized along the ridges the incoming electric field is perfectly reflected for a range of metal content. Conversely, in the IR-vis range the system is only partially reflecting having R ≈ 0.5, due to interband absorption. As for BMM mixture, the minimum of reflectance is at ∼ 2.2 eV for all metallic contents. Completely different is the optical response in the direction perpendicular to the ridges: R ⊥ is generally < 0.5 except for very high metal percentages (>80%). In particular, at low energies and low metallic content the material is essentially transparent, in agreement with the experimental findings 16 . Nevertheless, since the present theoretical description does not include the effect of the coupling with the lattice vibrations, we do not observe the Reststrahlen band and the modification of reflectance at 11µm, as discussed in the original paper 16 .
The different reflectivity of RHS is made evident from the comparison of the effective dielectric functions [panels (c),(d)] calculated for f = 0.5, which reproduces the experimental case of Ref. [ 16].˜ has the typical fingerprints of a d-metal: it is dominated by a Drude-like tail for E → 0 and by a crossover frequency of˜ r at E=0.75 eV, which corresponds to the excitation of a screened plasmon. Notably, with respect to the rutile phase, the plasma frequency is now strongly redshifted. This effect can be ascribed to insertions of the dielectric sections, which reduce the overall free electron density. The imaginary part,˜ i , is characterized by an absorption band in the IR-vis region due to many interband transitions, e.g. between occupied O(2p) or V (a 1g ) and the V (e π g ) states just above the Fermi Level, as pictorially indicated by red arrows in Fig. 2. Along the perpendicular direction, the dielectric function represents a semiconducting system, with˜ r converging to the static dielectric constant 0 = 25 for E → 0. Again, the presence of the metallic component increases the electrostatic screening to result in a higher value of 0 with respect to the M 1 phase and in the shrinking of the absorption band that now is singly peaked around E=0.75 eV, i.e. the energy value corresponding to the plasma frequency of the parallel component.
The directional optical response of this ridged structure indicates that ordered joint-phase heterostructures may sustain radiation with hyperbolic dispersion. In order to get deep insight on this point, we considered the standard geometry exploited in the fabrication of HMMs given by the alternate stacking of metallic and dielectric layers [ Fig. 4(c)]. Actually, this kind of VO 2 structure has not been fabricated, yet. However, the rapid and continuous improvements in the spatially selective control of metallic and semiconducting phases (e.g. checkerboard pattern, ridged and herringbone structures, stacked nanobeams) 52-55 make the growth of VO 2 layered heterostructures (LHS) very reasonable in the shortmedium term.
Assuming the stacking direction as the anisotropy optical axis, the parallel and perpendicular components of the effective dielectric function are given by 3 : where d m and d d are thicknesses of the metallic and dielectric layers [ Fig. 4(c)]. The RHS and the LHS model are formally equivalent as˜ (RHS) =˜ ⊥ (LHS) and We analyzed the hyperbolic dispersion condition of VO 2 LSHs, plotting the sign function Black areas in Fig. 7(a) identify the regions where the sign function S is negative, i.e. where the metamaterial has a hyperbolic behavior: two broad interconnected regions at 0-0.7 and 0.7-1.7 eV; one narrow region at 1.9-2.1 eV. To investigate the characteristics of these regions we select two values of metal content, namely f = 0.25 and f = 0.5 (labeled 1 and 2, respectively in Fig. 7) and we plot the corresponding parallel and perpendicular components of the dielectric function [panel (b)]. From the comparison between the two panels it follows that each area of panel (a) corresponds to different electrical characteristics. White sections, for which the dielectric functions have the same sign, correspond to an overall dielectric character, being both components of˜ r finite and positive. Black areas correspond to distinct hyperbolic types: at lower and higher energies (blue areas in Fig.  7) we see˜ > 0 and˜ ⊥ < 0, which is the definition of type-II HMMs, while in the intermediate region (orange areas) the condition˜ < 0 and˜ ⊥ > 0 identifies a type-I HMM. The energy alternation of these electrical behavior depends on the specific metal/semiconductor ratio in the sample. At low metal content (1) the hyperbolic regions are separated by dielectric energy gaps. In case 2, the sample switches from type-II to type-I at E=0.75 eV, where both components are zero (epsilon-near-zero ENZ regime 56 ).
Along the lines of Ref. [ 6], we define two parameters in order to quantify the hyperbolic character of LHS, namely the quality factor Q , and the strength of dielec-tric anisotropy ∆˜ : with j = , ⊥. The results for the key LHS systems 1 and 2 are summarized in Fig. 7(c). Although plotted over the entire energy range 0-3 eV for completeness, the quality factor Q accounts for the energy losses only along the direction in which the material has a metallic character. Thus, for the perpendicular component it is meaningful in the type-II region (blue area), while for the parallel component in the type-I region (orange area). Systems with hyperbolic character (S < 0) are considered good HMM if the maximum value of the quality factor is Q max > 3. Both systems satisfy this condition at low energy (E < 0.75 eV): they reach the maximum at E ∼0.3 eV, where Q max (1)=3.1 and Q max (2)=4.1. This qualifies VO 2 heterostructures as good type-II HMM materials in the THz and mid-IR regions. At higher energies, despite a formal hyperbolic character of the system, the metal energy loss prevents any realistic application of these LHSs as good metamaterials. The same conclusion comes from the analysis of the strength of dielectric anisotropy [ Fig. 7(d)]. Since the hyperbolic behavior derives from the spatial anisotropy of the optical response, the quality of HMM derives not only from the magnitudes of the constituent dielectric functions, but from their differences or ratios. As the metal ratio in the heterostructure increases, the maximum value of ∆˜ increases and diverges for E → 0 while its energy position is monotonically redshifted. This corresponds to the fact that for E → 0, the real part of the dielectric function in the metal component becomes so negative that there is an impedance mismatch with the other medium. Systems with ∆˜ max >20-30 are considered good HMMs 10 . The complete analysis of our results show that in the case of layered VO 2 , this condition is reached for f ≥ 0.25, with ∆˜ max ranging from mid-IR to THz regions. The strength of dielectric anisotropy also affects the dispersion relation of the radiation that may propagate along the medium, as and ⊥ define the geometric characteristics of the hyperbolic isosurface equation Let us consider the radiation pattern in the case of VO 2 LHSs: In general, electromagnetic waves cannot propagate in media with negative permittivity (e.g. metals). In the particular case of HMMs, the condition ⊥ < 0 implies that the system behaves like a metal in one direction and as a dielectric in the other. As a consequence, waves with arbitrarily large wavevectors k may have a propagating nature, while in isotropic materials they become evanescent due to the bound isofrequency contour (a sphere). In this sense, HMMs can be considered as plasmon-polaritonic crystals where the coupled states of light and electron density (i.e. plasmons) give rise to travelling extraordinary (TM) 57 waves, known as volume plasmon polaritons (VPP) 6,36 . To characterize these VPPs, we consider the electromagnetic field condition corresponding to the vector diagram shown in Fig.  8(a), where E is the electric-field vector, D is the electricdisplacement vector, and S is the Poynting vector. D lays in the principal plain containing both optical axis (aligned along z) and the wavevector. ϕ is the angle between the wavevector and the optical axis. In anisotropic materials, the vectors E and D are not usually parallel. An immediate consequence of this is that the Poynting vector S, which points in the direction of energy flow, and wavevector k, directed along the wavefront normal, need not to be parallel. The spatial anisotropy of HMM can be expressed through an angular dependence of the dielectric function 6 : The condition Re[ (ϕ c )] = 0 determines the angular boundary between the metallic and dielectric response of the metamaterial. The critical angle ϕ c is the direction of the incident radiation that may excite the formation of a travelling plasmon-polariton wave in the system. If we consider the interface with a another dielectric medium (e.g. air), the wavevector at the interface follows Snell's law. Thus, the angle θ between the extraordinary wave (i.e. the Poynting vector) and the optical axis is tan(θ) = Re[˜ ⊥ ]tan(ϕ) 6 . At the critical angle ϕ c , θ c reduces to the asymptotic expression θ c = arctan −Re[˜ ⊥ ] , and represents the angle of propagation of the VPP along the metamaterial. This means that the plasmon-polariton propagates throughout the volume of the HMM along a cone with axis coincident to the optical axis and angle θ c [ Fig. 8(a)]. The number of VPP modes that can be excited depends on the actual number of metal/dielectric layers in the heterostructure. The maximum number of modes in an ideal structure is always one less than the number of metallic layers. However, a real structure always has a maximum wavevector beyond which the VPP modes cannot exist. This happens when the size of the VPP wavevector becomes comparable to the size of the layers comprising the HMM and the modes no longer experience a homogeneous effective medium. The latter is also the validity breakdown of the effective medium approximation. The k wavevector follows a truly hyperbolic surface (i.e. open form) only for the ideal case with zero absorption processes. In the case of real absorbing materials, the isosurface equation reduces to a closed form, whose final dispersion crucially depends on the imaginary part of the dielectric function. The angular dependence of (ϕ) and the critical angle θ(ϕ) in the case of VO 2 LHSs (f =0.5) are shown in Fig. 8(b-c), for different values of the incoming radiation energy. In the low energy region (E < 0.75 eV), where the system has a type-II character (black and red lines), the absolute value of (ϕ) decreases, while ϕ c shifts to larger angles (ϕ c = 28 • and 57 • ) as the value of E is increased [panel (b)]. The dielectric functions look like a single dipolar resonance with a high value of the imaginary components (not shown). The corresponding negative critical angles [θ c = −61 • and −31 • , panel (c)] indicate the propagation of lateral backward wave with respect to the interface without negative refraction 57 . At E = 0.75 eV the real part of (ϕ) is zero and θ c is indefinite. Approaching this critical value in the ENZ region, e.g. E ∼ = 0.75 eV (blue line), (ϕ) crosses the zero at ϕ c = 68 • . The corresponding critical angle θ c is very small, imparting an extreme conical focusing to the propagation wave. At E = 1.0 eV the system has a type-I character (green line) with a small valued dielectric function, typical of dielectric-like HMMs. In this case the critical angle θ c is small and positive, which corresponds to a standard positive wave refraction within a collimated propagation cone, and could be exploited for applications in hyperlenses.

IV. CONCLUSION
We presented a first principles investigation of the optical properties of VO 2 and used EMA to determine the properties of mixed-phase metamaterials obtained by combining different phases of the same VO 2 compound, i.e. without the introduction of other different media, to form structured materials with original and tunable optoelectronic properties and (in certain cases) with highly anisotropic permittivity. Our results complement and agree with experiments performed on both disordered mixtures and ordered ridged structures. Furthermore, we demonstrated the possibility to use a variable amount of metallic VO 2 content in device applications where it is necessary to control the optical properties at a target wavelength.
Finally, we investigated the hyperbolic behavior and the plasmonic characteristics of stacked VO 2 heterostructures, a still unexplored possibility made feasible by the unique properties of this phase change compound. The possibility to realize joint-phase mixtures allows one to combine the optical constants of single phases to form a ready-made metamaterial with unforeseen properties. Indeed, this combines the structural tunability typical of artificial HMMs with reduced intermaterial interfaces typical of natural HMMs, and thus opens the way for a new class of metamaterials with controllable optical properties over a wide energy range (THz to visible).

V. ACKNOWLEDGMENTS
We thank Thushari Jayasekera for her friendly support and advice. ME acknowledges the Southern Illinois University Study Abroad Program for facilitating international collaboration in addition to the SIU Chancellor's Scholar Program and Center for Undergraduate Research and Creative Activities (CURCA) for partial financial support.