Fuyang Tay

Ali Mojibpour

Stephen Sanders

Shuang Liang

Hongjing Xu

Geoff C Gardner

Andrey Baydin

Michael J Manfra

Alessandro Alabastri

David Hagenmüller

Junichiro Kono
email: kono@rice.edu.

Ultrastrong photon-photon coupling

Recent studies have shown that matter can ultrastrongly couple with the quantum vacuum field inside a photonic cavity, producing a nonclassical ground state that contains a finite number of photons. Here, we present a novel mattervacuum hybrid in a multimode photonic cavity whose ground state contains ultrastrongly coupled photons. This unique photon-photon coupling was realized in a three-dimensional terahertz photonic-crystal cavity, where two adjacent cavity modes mixed together through simultaneous coupling with the cyclotron resonance of a two-dimensional electron gas with a coupling strength exceeding the intermode frequency. Our microscopic theory successfully explains the salient features of our experimental observations, highlighting the spatial overlap of mode profiles as a key enabler of photon-photon ultrastrong coupling. Our findings provide guidelines for harnessing photon-photon correlations for furthering the physics of vacuum-dressed matter as well as for developing vacuum-enabled quantum technology.

Introduction

Correlations between photons play a crucial role in fundamental quantum optics research and quantum technology, such as nonclassical light generation and quantum information processing (1,2). Generating such a state, including squeezed states (3)(4)(5), requires external driving, and hence, the existence of such a state is only transient. Recent theoretical studies have suggested that significant correlations between photons and matter excitations can exist in a cavity quantum electrodynamics (cQED) system without external driving. Specifically, when the light-matter coupling strength, g, is comparable to the bare frequencies of light and matter, the ground state is predicted to be a two-mode squeezed vacuum (6,7). This regime, known as the ultrastrong coupling (USC) regime of light-matter interaction (8)(9)(10), has been achieved in several solid-state platforms, paving the way to realizing macroscopic ground-state squeezing.

Here, we introduce the novel concept of ultrastrong photon-photon coupling in a multimode cQED system. Previous studies have shown that, when g becomes comparable to, or larger than, the spacing between adjacent cavity modes, ∆ω, in a multimode cavity, the system enters the superstrong coupling (SSC) regime (also known as the multimode strong-coupling regime (11)), where matter exchanges energy with the cavity faster than the photon roundtrip time (12,13).

In this regime, intermode coupling is expected, but the ground-state properties of a system that simultaneously achieves USC and SSC have not been explored. An interesting question is whether USC between different photonic modes i and j can be achieved, which should lead to finite intermode correlations in the ground state, i.e., ⟨G|â † i âj |G⟩ > 0, where â † and â are the photon creation and annihilation operators and |G⟩ represents the ground state. Such multimode ground-state correlations would open new possibilities for harnessing vacuum electromagnetic fields to generate multimode nonclassical states of light (14,15) and for studying the many-body regime of cQED in the presence of nonlinearities (16)(17)(18).

We realized USC between two photonic modes of a THz three-dimensional photonic-crystal cavity (3D-PCC). The coupling was mediated by the cyclotron resonance (CR) of a Landauquantized two-dimensional electron gas (2DEG) in GaAs quantum wells (QWs). Namely, the 2DEG CR ultrastrongly coupled with both photonic modes, achieving USC and SSC simultaneously. The cavity had an asymmetric design such that the mode spatial profiles in the 2DEG region depended on the THz electric-field polarization. The magnetic-field dependence of the experimentally observed Landau polariton frequencies was well described by a microscopic quantum model following Ref. (19). Despite the lack of nonlinearity in the GaAs 2DEG, the calculations suggest that the intermode correlations in the ground state, ⟨G|â † 1 â2 |G⟩, was comparable to the intramode correlations (or numbers of photons) in the ground state ⟨G|â † 1 â1 |G⟩ and ⟨G|â † 2 â2 |G⟩. We show that the full in-plane spatial variation of the cavity electromagnetic field has to be considered to reproduce the observed polariton spectra. Although g > ∆ω was always satisfied, intermode correlations in the ground state occurred only when the mode profiles significantly overlapped spatially. We propose a relevant figure of merit (FOM) to quantify photon-photon couplings in multimode systems, which governs the correlations between photonic modes when the matter and photon frequencies are close. We found the conditions for achieving and controlling ultrastrong photon-photon coupling, which show that the 3D-PCC is an ideal candidate for exploring the many-body regime of cQED in the THz frequency range.

Defect-engineered three-dimensional photonic-crystal cavity

Let us consider a cQED system where N cavity modes interact with a single matter excitation.

When g < ∆ω, vacuum Rabi splitting (VRS) occurs at N well-separated values of the detuning, as schematically shown in Fig. 1A for N = 2. When the system enters the SSC regime, g > ∆ω, the matter mode can couple with multiple cavity modes simultaneously at a single value of the detuning. However, if the entire cavity is filled with matter, different photonic modes separately couple with the matter at different spatial locations because the electric-field profiles of the cavity modes are orthogonal to each other over the entire volume of the cavity (assuming no dissipation) (20,21). Consequently, no matter-mediated intermode coupling exists, and each photonic mode splits into an upper polariton (UP) and a lower polariton (LP), resulting in a total of 2N branches. This situation, depicted in Fig. 1B for N = 2, is referred to as the 2N scenario (21,22).

To induce coupling between cavity modes, one can reduce the active region, i.e., the volume inside the cavity occupied by the matter (20,21). This is because the spatial profiles of different cavity modes over a portion of the volume inside the cavity are typically not orthogonal. When the spatial profiles of different cavity modes within the active region highly overlap, these modes can couple with each other through the matter, producing N + 1 polariton branches with an Sshaped middle polariton (MP) (20); see Fig. 1C for N = 2. This situation is known as the N + 1 scenario.

We used a 3D photonic crystal, i.e., a dielectric structure that has periodic modulations of refractive index in all three dimensions and exhibits a complete photonic band gap. Adding a defect into the 3D photonic crystal breaks discrete translational symmetry, inducing photonic modes within the band gap (23). The electromagnetic fields at these photonic mode frequencies are confined in the vicinity of the defect, which can be engineered to tailor the mode spatial profiles. In addition, multiple cavity resonances exist closely in frequency within the photonic band gap, and hence, a THz 3D-PCC provides a suitable platform for exploring multimode USC through coupling with a collective matter excitation confined in the defect.

We embedded a multiple-QW 2DEG layer (∼2 µm thick) on a 58-µm-thick planar GaAs substrate between a pair of woodpile structures, often adopted 3D photonic crystals created by stacking rod arrays periodically with alternating orientations (Fig. 1D). There is no translational invariance along the stacking direction, z, while (discrete) translational symmetry exists in the x-y plane. We computed the projected band structure of the bare woodpile cavity without the 2DEG layer, with infinite periods above and below the GaAs defect, in the in-plane irreducible first Brillouin zone using plane-wave expansion (Fig. 1E). The dimensions of the woodpile structure are shown in Fig. 2A. Guided modes (black lines in Fig. 1E) that are exponentially localized in the vicinity of the defect (along the z direction) and propagate in the x-y plane are present in the photonic band gap (between 318 GHz and 478 GHz at the Γ point) (24).

Our designed woodpile structure contained only two periods on each side of the cavity in the z direction, so the modes in the photonic band gap had finite lifetimes. The lifetimes were estimated from transmission spectra via numerical simulations, as shown in Fig. 1F, where the 3D-PCC was probed at normal incidence (corresponding to the Γ point). Our 3D-PCC hosted modes with varying quality factors depending on the position of the modes within the band gap: modes with frequencies, ω 1 /2π = 338 GHz, ω 2 /2π = 382 GHz, ω 3 /2π = 417 GHz, and ω 4 /2π = 468 GHz exhibited quality factors of 72, 70, 6200, and 1540, respectively. We used deep reactive ion etching to create a bar array in a silicon wafer (Fig. 1G). The final woodpile structure was formed by mechanically stacking multiple patterned silicon wafers.

The active region for light-matter coupling in our system is restricted to the 2D plane occupied by the 2DEG. Our cavity structure lacks mirror symmetry in the vertical direction because the silicon bars at the edges of the defect are oriented perpendicular to each other (Fig. 2A).

When there is no 2DEG layer inside the cavity, the transmittance spectrum is the same for xand y-polarized THz light (24,25). However, even though they have the same frequencies, their in-plane mode profiles at the vertical location of the 2DEG, z 2DEG , vary, depending on the polarization of the incident THz radiation. We refer to the modes excited by x-and y-polarized incident THz radiation as the transverse-magnetic (TM) and transverse-electric (TE) modes, respectively. The electric energy density of cavity mode 1 is maximized at slightly different vertical locations for the TM and TE modes (Fig. 2, B and C) (24). Figures 2, D andE, display the spatial profiles of modes 1 and 2 at z 2DEG for the TM modes, in which the cavity electric field is tightly confined at the edges of the bar. By contrast, the mode profiles for the TE modes distribute more evenly within the unit cell (Fig. 2, F andG). Therefore, the spatial overlap of the mode profiles within the active region is different for the TE and TM modes.

Polarization-dependent mixing of photonic modes

Figures 3, A andB, show transmittance spectra for our 3D-PCC with a 2DEG layer as a function of magnetic field applied along the z direction, B, for the TM and TE modes, respectively.

The white dots denote the B-dependent polariton frequencies that were extracted from longer time-domain traces (24). For simplicity, we restrict ourselves to the first two modes p = 1, 2 in the following discussion. Three polariton branches, one UP and two LPs, are observed for TM (Fig. 3A), accompanied by a splitting spanning the space around the line of the bare CR frequency (white dashed), ω c , that separates the UP and LPs. This is consistent with the 2N scenario (Fig. 1B). By contrast, such a splitting is not visible for TE (Fig. 3B). Instead, a polariton branch that crosses the CR line and resembles the S-shaped MP in the N + 1 scenario (20,22) (Fig. 1C) is observed.

To gain a better understanding of the observed polarization-dependent transmittance spectra, we developed a microscopic quantum model described by the Hamiltonian Ĥ = Ĥcav + ĤCR + Ĥint + ĤA 2 (24). The free-photon Hamiltonian is Ĥcav = p,j ℏω p â † p,j âp,j , with âp,j the annihilation operator of a photon in mode p with polarization j and frequency ω p . For simplicity, we did not introduce another index in the photon operators to distinguish between the

The next term in the Hamiltonian

Ĥint = p,j dρ a ℏg p,j (ρ) b(ρ) + η j b † (ρ) âp,j + η j â † p,j (1)
is the linear coupling between the CR and each cavity mode (p, j), whose strength is g p,j (ρ) = E p,j (ρ, z 2DEG ) e 2 ω c n e /(4ε 0 m eff ω p a). Here, η x = +, η y = -for the two cavity photon polarizations, e and m eff are the electron charge and mass in GaAs, respectively, ε 0 is the vacuum permittivity, and n e is the total electron density. Since the electric field in the 3D-PCC is not uniform in the plane (in contrast to 1D-PCCs (26,27)), the coupling of CR excitations to each mode (p, j) is weighted by the (real, dimensionless) electric field mode functions E p,j at the in-plane position ρ and the vertical location z 2DEG of the QWs.

The final term in the Hamiltonian, the so-called A 2 term

ĤA 2 = p,p ′ ,j dρ a 2 ℏg p,j (ρ)g p ′ ,j (ρ)
ω c âp,j + η j â † p,j â † p ′ ,j + η j âp ′ ,j (2)
is responsible for direct intermode coupling. This term scales with dρ g p,j (ρ)g p ′ ,j (ρ)/ω c , and thus, involves an overlap integral between the in-plane spatial profiles of the modes E p,j (ρ, z 2DEG ).

As mentioned previously, the cavity modes over the entire volume inside the cavity are orthogonal to each other, dρdz ε(ρ, z)E p,j (ρ, z)E p ′ ,j ′ (ρ, z) = a 3 δ p,p ′ , with ε(ρ, z) the inhomogeneous dielectric profile. However, the cavity mode profiles within the 2D plane at z 2DEG do not necessarily follow the orthogonality relation (20). Alternatively, the orthogonal modes can also be quantum mechanically coupled through dissipation, which can be described by quantized quasinormal mode theory (28).

The lossless, multimode Hamiltonian, Ĥ, can be simplified to a decoupled Hamiltonian under certain circumstances. We introduce new CR excitation operators following the spatial profiles of the cavity modes, bp,j = (dρ/a) (g p,j (ρ)/ g p,j ) b(ρ), with the effective coupling strength g p,j = (dρ/a 2 ) g 2 p,j (ρ). The commutation relations between these new CR modes can be put in the form [ bp,j , b † p ′ ,j ′ ] = ξ p,p ′ ;j,j ′ , where ξ p,p ′ ;j,j ′ = (dρ/a 2 ) g p,j (ρ)g p ′ ,j ′ (ρ)/( g p,j g p ′ ,j ′ ) (0 ≤ ξ p,p ′ ;j,j ′ ≤ 1) is proportional to the in-plane spatial overlap of the different cavity modes, exactly as the off-diagonal contributions to the A 2 term, and is similar to the normalized overlap factor defined in Ref. (20). Assuming that the new CR modes are orthogonal, which is a priori not the case, one would have ξ p,p ′ ;j,j ′ = δ p,p ′ δ j,j ′ . The full Hamiltonian would then be simplified to a "decoupled" Hamiltonian with ĤCR = ℏω c b † p,j bp,j , and

Ĥint = ℏ g p,j bp,j + η j b † p,j âp,j + η j â † p,j ĤA 2 = ℏ g 2 p,j ω c âp,j + η j â † p,j â † p,j + η j âp,j . (3)
The deviation of the parameter ξ p,p ′ ;j,j ′ with respect to δ p,p ′ δ j,j ′ , therefore, tells us how accurate the approximation of using the decoupled Hamiltonian to compute the polariton branches is as compared to the full Hamiltonian. For instance, ξ p,p ′ ;j,j = 1 (ξ p,p ′ ;j,j = 0) corresponds to perfect overlap (no overlap) between the cavity modes p and p ′ both with polarization j, in which case the cavity modes are coupled (decoupled) via the CR.

We computed transmittance spectra by extending the input-output model of Ref. (29) in a simple planar geometry to a multimode PCC cavity (24). We took into account dissipation through the (Markovian) coupling of cavity modes and CR excitations to phenomenological bosonic reservoirs, with associated quality factors of the bare cavity modes and the intrinsic CR decay rate, respectively (24). For our 3D-PCC system, there is only a small overlap between the spatial profiles of the TM modes 1 and 2 (ξ 1,2;x,x = 0.29), suggesting that the intermode coupling is weak. The decoupled model in Eq. ( 3) is thus expected to provide a faithful description of the system and in particular the existence of a LP-UP splitting close to resonance, which is consistent with the transmittance spectra calculated using the full and decoupled Hamiltonian in Fig. 3, C and E. Conversely, for TE modes, the spatial profiles of cavity modes 1 and 2 highly overlap (ξ 1,2;y,y = 0.91), suggesting that they are strongly coupled through the CR.

The full model predicts the emergence of an S-shaped MP that adiabatically transforms from the UP of mode 1 to the LP of mode 2, passing continuously through the CR line (Fig. 3D).

Note that the small LP-UP polariton splitting would vanish if ξ 1,2;y,y is exactly equal to 1 (24).

The decoupled model that predicts the LP-UP splitting thus exhibits large deviations from the experimental data points (Fig. 3F). In addition, we find that ξ p,p ′ ;x,y < 2 × 10 -3 for all p, p ′ , indicating that cavity modes with different polarizations are weakly coupled.

Ultrastrong photon-photon coupling

The strength of the light-matter interaction is usually estimated using the on-resonance VRS obtained from experiments, which is typically equal to twice the coupling constant of the linear coupling term Ĥint . However, in the 3D-PCC, the VRS at zero detuning cannot accurately determine the strength of the light-matter interaction for two reasons: (i) the coupling constant depends on the position ρ in the unit cell, while the VRS is contributed by the coupling at all ρ;

and (ii) the shift of the polariton branches, which is induced by the intermode coupling, reduces the VRS (e.g., the redshift of the UP marked by the yellow arrows in Fig. 3, D andF).

In order to quantify the strength of all possible couplings for a multimode system in the USC regime, we propose to use a generalized FOM defined as

η pp ′ ,j ≡ (dρ/a 2 ) g p,j (ρ)g p ′ ,j (ρ) ω c (ω p + ω p ′ )/2 . ( 4
)
Note that η pp ′ ,j does not depend on B, and thus it is independent of detuning. The strength of the coupling between the CR and the cavity mode p can be characterized by the "diagonal" contributions (p = p ′ ), which is the same as the conventional FOM for the USC regime (8-10), except that it includes the integral of g p,j (ρ) over the whole unit cell. For TM modes, η 11,x ≈ η 22,x = 0.21, η 33,x = 0.22, and η 44,x = 0.19, while for TE modes, η 11,y = 0.2, η 22,y = 0.17, η 33,y = 0.08, and η 44,y = 0.19. Hence, our system reached the USC regime for both the TE and TM modes since η pp,j ≳ 0.1 for all p and j. The "off-diagonal" contributions (p ̸ = p ′ ), which are related to the spatial overlap of the mode profiles (ξ p,p ′ ;j,j ∝ η 2 pp ′ ,j ), determine the effective coupling between different cavity modes with the same polarization j. The parameter η pp ′ ,j serves as a good FOM to quantify the photon-photon coupling because the intermode ground-state correlations scale with it (24). As η 12,x = 0.11 and η 12,y = 0.17, our system was in the USC regime in terms of photon-photon coupling (η pp ′ ,j ≳ 0.1).

We computed the photon ground-state (virtual) correlations ⟨â † p,j âp ′ ,j ⟩ for both the TE and TM modes by inverting the Bogoliubov transformation that diagonalizes the full (quadratic) Hamiltonian; see Fig. 4. The intermode correlation in the ground state for the TE modes, ⟨â † 1,y â2,y ⟩, lies between the intramode correlations in the ground state, ⟨â † 1,y â1,y ⟩ and ⟨â † 2,y â2,y ⟩, confirming the existence of significant correlations between the photonic modes in the ground state of the system. Furthermore, the intermode correlations for the TE modes (Fig. 4B) are larger than those for the TM modes (Fig. 4A) by a factor ≃ 2.6. Note that those intermode ground-state correlations are contributed by both Ĥint and ĤA 2 terms (24). Further, we used the standard FOM g p,j /∆ω pp ′ , with ∆ω pp ′ = |ω pω p ′ |, to demonstrate that our system was in the SSC regime. At resonance between the first cavity mode and the CR (ω c = ω 1 ), g 1,x /∆ω 21 = 1.71 and g 2,x /∆ω 21 = 1.84 for TM modes, while g 1,y /∆ω 21 = 1.61 and g 2,y /∆ω 21 = 1.46 for TE modes. Thus, g p,j /∆ω pp ′ > 1 for both TE and TM modes.

Discussion

We reported ultrastrong photon-photon coupling in a THz cQED system that simultaneously achieves USC and SSC, namely, a Landau-quantized 2DEG in GaAs coupled to multiple photonic modes of a THz 3D-PCC. We introduced a detuning-independent parameter, η pp ′ ,j , which can be generalized to a larger matrix for more modes, to quantify couplings in multimode systems. Our findings showed that the ground-state correlations between the photonic modes are controlled by the frequencies of both matter and photonic modes, as well as the spatial overlap of the cavity mode profiles. Due to the absence of mirror symmetry in the stacking direction of our cavity, the TE and TM photonic modes exhibited different in-plane mode profiles at the 2DEG position, leading to a switchable photon-photon coupling while preserving the USC between the matter and each photonic mode. Our results highlight the importance of considering the full spatial dependence of the electric field in modeling the light-matter interaction inside cavities with spatially nonuniform electromagnetic fields, even if the spatial variation of the mode profiles is almost negligible within the typical wavelengths associated with the coupled electronic degrees of freedom.

Although the matter excitation in the present work is linear in the sense that the energy spacings between adjacent Landau levels are constant, further studies on nonlinear systems can be conducted by replacing GaAs with a non-parabolic semiconductor or graphene. The ultrastrong photon-photon coupling opens up interesting perspectives to realize vacuum-induced effects (30)(31)(32)(33) and explore many-body regimes of quantum optics with quantum vacuum fields from different modes. The 3D-PCC design is highly versatile and can be tailored to achieve photonic modes with smaller mode volumes and higher quality factors (34). Our approach can be utilized to investigate the coupling between the topological edge states of a 2DEG and the topological edge states of a chiral woodpile structure (35,36). Moreover, the photonic modes with well-defined in-plane wavevectors in the 3D-PCC satisfy the condition to circumvent the no-go theorem for the Dicke superradiant phase transition (or photon condensation) (37)(38)(39).

While the photon wavevector ∼ 2π/a is limited by the lattice parameter a of the 3D-PCC and remains relatively small compared to the electron momentum ∼ 1/l c , the possibility of fulfilling the criterion for 2D photon condensation in this regime has been recently discussed (38).

S1.2 Fabrication of the THz woodpile cavity

A 100-µm-thick silicon wafer was coated by lift-off resist and photoresist layers through spin coating. The layers were patterned by photolithography. An ∼ 1.2-µm-thick Al 2 O 3 layer was deposited on the coated silicon wafer by an e-beam evaporator. The lift-off process was completed by removing the resist with Remover PG; leaving Al 2 O 3 on the patterned area. Afterwards, the silicon wafer was etched by deep reactive ion etching using a Bosch process. During this process, the silicon wafer was cut into multiple smaller pieces, and the gaps between silicon rods and holes at the corners were also created in each piece of the smaller wafer. The remaining Al 2 O 3 layer was removed by immersing the sample in a heated 1:3 solution of phosphoric acid and sulphuric acid. The width of the silicon rods and the periodicity were 87 µm and 333 µm, respectively. A custom-made sample holder was used for stacking the silicon wafers to form a woodpile cavity. The orientation of the silicon rods was fixed by the screws that passed through the holes at the corners. The woodpile cavity can be conveniently disassembled and reassembled.

S1.3 THz time-domain magnetospectroscopy measurements

Transmittance spectra for samples were measured by a home-built THz time-domain magnetospectroscopy setup (26,27,(40)(41)(42)(43)(44)(45)(46)(47)(48).

S2.2 Microscopic model

The microscopic Hamiltonian is derived in the Coulomb gauge by extending the models of Refs. (19,27,49). We use the vector potential operator

Â(ρ, z) = p,j ℏ 2ε 0 ω p a 3 E p,j (ρ, z)e j âp,j + â † p,j . (S1)
At normal incidence (Γ point of the first Brillouin zone), the electric field mode functions E p,j (ρ, z) can be chosen real without loss of generality. The calculation can be trivially extended to the full Brillouin zone using complex fields with a quasi-momentum q as computed by FDTD. The linear coupling term between light and matter reads

e m eff dρdz Ψ † (ρ, z) π • Â(ρ, z) Ψ(ρ, z), ( S2
)
with π the gauge-invariant in-plane momentum including the contribution of the static magnetic field B, and the fermion field Ψ(ρ, z) = L -1/2 n,k e -iky χ n,k (x)ξ(z)ĉ n,k (49). We use the Landau gauge, where single electron states are characterized by an integer n and the y component of the electron momentum k. The wave functions are plane waves in the y direction (L is the length of the 2DEG in the plane), and those of a harmonic oscillator (Hermite polynomial of order n) centered at the guiding center position k l 2 c , with l c = ℏ/eB the magnetic length. When computing the matrix elements entering Eq. (S2), we assume that E p,j (ρ, z) remains constant over the QW thickness, and exploit the discrete translational invariance in the plane by decomposing the electric field spatial functions into Fourier series E p,j (ρ, z 2DEG ) = G U p,j (G)e iG•ρ , with G = (2πm x /a)e x + (2πm y /a)e y (m x , m y ∈ N) the reciprocal lattice vectors. One is then left with the calculation of overlap integrals of the kind

I n ′ ±1,k ′ n,k (G) = dy L e i(k-k ′ +Gy)y dx χ n,k (x)χ n ′ ±1,k ′ (x)e iGxx . ( S3
)
While the integral in the y direction simply provides the selection rule k ′ = k + G y , the calculation of the integral along x is performed using the dipole approximation. The mode functions E p,j (ρ, z 2DEG ) are slowly varying over the typical extent of the harmonic oscillator wave functions χ n,k (x). For a magnetic field B ≃ 1 T, the latter is indeed ∼ l c ≃ 10 -100 nm, which is much smaller than the woodpile lattice parameter a = 333 µm. Since the Fourier coefficients U p,j (G) exhibit large peaks at G j ∼ 1/a (m x ∼ m y ∼ 1) and rapidly decrease as

G j a → ∞, one has G j l c ∼ l c /a ≪ 1. Equation (S3) thus provides I n ′ ±1,k ′ n,k (G) ≈ δ k ′ ,k+Gy δ n,n ′ ±1 e iGxkl 2 c .
We introduce the CR excitation creation operator The transmission spectra shown in the main text were computed using an input-output model, which is an extension of the one introduced in Ref. (29) in a simple planar geometry.

b † (ρ) = 1 a √ N k,G ĉ † ν,k-Gy ĉν-1,k e
Transmission of THz radiation through the 3D-PCC is modeled by introducing two identical photon reservoirs on each side of the cavity in the z direction (top, bottom). Similarly, CR excitations of the 2DEG acquire a finite lifetime by interacting with a phenomenological bosonic reservoir. The total Hamiltonian includes the contribution (neglecting counter-rotating terms) ĤR = p,j,λ dq ω p,j (q)α † p,j,λ (q)α p,j,λ (q) + iκ p,j (q) αp,j,λ (q)â † p,jâp,j α † p,j,λ (q)

+ dq dρ ω(q, ρ) β † (ρ, q) β(ρ, q) + i κ(ρ, q) β(ρ, q) b † (ρ) -b(ρ) β † (ρ, q) . (S4)
The first contributions in the first and second lines of Eq. (S4) describe the energy of the photonic modes in each reservoir λ = (top, bot) and that of the matter excitation reservoir, respectively. q denotes a phenomenological (continuous) parameter ensuring that the reservoirs have continuous spectra. The second terms in both lines describe the coupling between the mode of the reservoirs and the modes of the system, which provides the latter with a finite lifetime. The input-output method consists in solving the equations of motion of the system by introducing "input" and "output" operators: αin p,j,λ (q) = lim t→-∞ αp,j,λ (q, t)e iω p,j (q)t , αout p,j,λ (q) = lim t→+∞ αp,j,λ (q, t)e iω p,j (q)t , and similar expressions for the matter excitation reservoir operators. One can then find a linear relation between the input and output operators, which allows to compute the normalized transmission spectrum (transmitted light with polarization j ′ , incident light with polarization j) as

T j,j ′ (ω) = p ′ ⟨α out † p ′ ,j ′ ,bot (q)α out p ′ ,j ′ ,bot (q)
⟩ p ⟨α in † p,j,top (q)α in p,j,top (q)⟩ .

We assume THz radiation coming from the top reservoir and populating all photonic modes equally, i.e., ⟨α in † p,j,top (q)α in p,j,top (q)⟩ does not depend on p nor j, while ⟨α in † p,j,bot (q)α in p,j,bot (q)⟩ = 0 ∀ (p, j) and ⟨ βin † (ρ, q) βin (ρ, q)⟩ = 0 ∀ ρ. The output expectation value ⟨α out † p,j,bot (q)α out p,j,bot (q)⟩ depends on the decay rates of each photonic modes Γ p,j = ω p /Q p,j ≡ πκ 2 p,j (q 0 )ρ(ω), with Q p,j the quality factor of the photonic mode (p, j) as computed in FDTD (MEEP implementation ( 50)) without the 2DEG, and on the CR excitation decay rate Γ c ≡ π κ 2 (ρ, q 0 )ρ(ω). The latter corresponds to the intrinsic CR decay rate. The quality factors of the photonic modes were obtained by fitting the peaks in transmittance spectra with Lorentzian functions. However, deviations from pure Lorentzian line shapes were clearly observed (especially for low-Q modes), similar to a previous report for 0D-PCCs (51), which has been attributed to dissipation (radiative leakage of the cavity modes) [START_REF] Sauvan | Theory of the spontaneous optical emission of nanosize photonic and plasmon resonators[END_REF]. Note that the CR decay rate is suppressed when the 2DEG is placed inside a high-Q cavity (26). As a good estimate, we used the CR decay rate obtained from previous experiments with a 1D-PCC (Γ c /2π = 5.7 GHz) (26). In the previous equations, q 0 and q 0 are solutions of ωω p,j (q) = 0 and ωω(q, ρ) = 0, respectively, and ρ(ω) is the effective density of states of the reservoirs, which are assumed to be Markovian, i.e., ρ(ω) is supposed to be frequency-independent in the range of interest.

In order to further characterize the full and decoupled models discussed in the main text, we introduce a toy model with two photonic modes, only one polarization, and continuously tunable spatial overlap between the photonic modes. The toy model Hamiltonian reads

H = p ℏω p â † p âp + dρℏω c b † (ρ) b(ρ) + p dρ a ℏg p (ρ) b(ρ) + b † (ρ) âp + â † p + p,p ′ dρ a 2 ℏg p (ρ)g p ′ (ρ) ω c âp + â † p â † p ′ + âp ′ , (S5)
with the coupling strength g p (ρ) = E p (ρ) e 2 ω c n e /(4ε 0 m eff ω p a), and the in-plane mode pro-files

E 1 (ρ) = sin 2πx a sin 2πy a E 2 (ρ) = sin 2πx a + (1 -ϵ)π 2 sin 2πy a + (1 -ϵ)π 2 . ( S6

S2.3 Ground-state correlations from different terms in Hamiltonian

Figure S3 shows the contribution of the ground-state correlations from different terms in the Hamiltonian. The Ĥint term provides the magnetic field dependence of the ground-state correlation (Fig. S3A). The ĤA 2 term provides the baseline of the ground-state correlation with the full Hamiltonian (Fig. S3B), supporting the claim that the η pp ′ ,j is a suitable FOM for ultrastrong photon-photon coupling. The ground-state correlations are zero if the antiresonant part of the Ĥint and ĤA 2 terms are neglected (Fig. S3C). The ground-state correlation will be overestimated if the resonant part of the Ĥint and ĤA 2 terms are neglected (Fig. S3D).

S2.4 Vacuum fluctuations and localization of the cavity modes

The vacuum fluctuations of the photonic modes shown in the main text were computed as follows: The variance of the electric energy density along z in the vacuum state |0⟩ reads

I(z) = dρ ε 0 ε(ρ, z) ⟨0| Ê2 (ρ, z) |0⟩ /2,
E p (ρ) = ℏω p 2ε 0 a 3 × E 2 p,x (ρ, z 2DEG ) + E 2 p,y (ρ, z 2DEG ). (S9)
In order to assess the degree of localization of the different cavity modes in the z direction, one can treat the electric energy density of mode (p, j) as a normalized probability distribution, i.e., (about two Si logs thick), and σ 3,x ≈ σ 4,x ≈ σ 3,y ≈ σ 4,y = 0.2a (less than one Si log thick). For all cavity modes, more than 97% of the electric energy is located in the 3D-PCC. As expected, the localization of the cavity modes along the z direction increases as the latter are located deeper into the photonic band gap.

S2.5 Numerical simulations

The mode profiles of the cavity modes induced by a 60 µm-thick bare GaAs layer were calculated by the MEEP software (50). The electric field was assumed to be constant over the thickness d ≃ 2 µm of the QW heterostructure. The transmittance spectrum of the bare cavity in the main text was simulated by using the COMSOL Multiphysics software. We used the permittivity ε = 11.6964 for the silicon layers, and ε = 12.96 for the GaAs layer. For the simulations with the 2DEG layer (Fig. S4), because the woodpile structure exhibits mirror symmetries, only 1/4 of the unit cell of the woodpile structure was considered in the geometry to reduce the simulation time. A transition boundary condition with an effective thickness, d = 2 µm, was used to emulate the MQW structure. A gyrotropic permittivity tensor was used to describe the complex permittivity of the 2DEG layer at different magnetic fields:

ε =   ε xx (ω) ε xy (ω) 0 -ε xy (ω) ε xx (ω) 0 0 0 ε zz   , (S13)
with

ε xx (ω) = ε bg - ω 2 pl (ω -iγ) ωd[(ω -iγ) 2 -ω 2 c ] , (S14)
ε xy (ω) = -iω 2 pl ω c ωd[(ω -iγ) 2 -ω 2 c ] , (S15)
ε zz = ε bg , (S16)
where ε bg = 12.96 is the background permittivity, ω pl = n e e 2 /(ε 0 m eff ) is the plasma frequency and ω c = eB/m eff is the cyclotron frequency.

S2.6 Thickness of the defect layer in numerical calculations

Due to imperfections in mechanical polishing, the thickness of the GaAs substrate of the 2DEG layer in the experiment was not perfectly homogeneous. To determine the suitable thickness of the substrate for numerical simulations, we ran simulations in the presence of the 2DEG layer at B = 0, 1.5, 2.5, and 7 T with different thicknesses of the GaAs substrate. We calculated the deviations of the polariton frequencies for TE and TM modes between the experiment and simulation, ∆f P,B = f expf sim , where P = {UP, LP1, LP2} at each B. Finally, we calculated the root mean square of the total deviation,

D = TE,TM B P ∆f 2 P,B 2N , (S17)
where N = 13 is the number of data points. The optimized thickness of the defect layer that gives the minimum deviation is 60 µm (including the 2-µm-thick MQW layer), as shown in Fig. S5. This is consistent with the profilometer measurement.

S2.7 Simulations with different ranges of time delays

We investigated the transmittance spectra of the system with different ranges of time delays in FDTD simulations (Lumerical). As shown in Fig. S6, in order to resolve modes 3 and 4 with high Q-factors in the spectrum, a long time-domain trace is needed. The amplitudes of modes 3 and 4 are weak for a time-domain trace up to 200 ps.

S2.8 Extraction of peak frequencies

The linewidth of the polariton branches in the color plot in the main text is limited by the frequency resolution of the measurements as a short time-domain range (33 ps) was considered.

Scanning a long range of time delays is required for THz-TDS measurements to resolve the actual linewidth of the coupled modes in the high-Q 3D-PCC. For instance, Fig. S7,A Echoes of THz pulses come from multiple optical components in the system, for instance, cryostat windows and crystals.

  TE and TM modes; instead, we directly mention it in the text. The effective CR Hamiltonian ĤCR = ℏω c b † (ρ) b(ρ) is written in terms of the collective excitation operators between adjacent Landau levels b(ρ) and b † (ρ) at the in-plane position ρ, which approximately satisfy the bosonic commutation relations [ b(ρ), b † (ρ ′ )] = δ(ρρ ′ ).

  52. C. Sauvan, J. P. Hugonin, I. S. Maksymov, P. Lalanne, Theory of the spontaneous optical emission of nanosize photonic and plasmon resonators, Phys. Rev. Lett. 110, 237401 (2013).

Fig. 1 :Fig. 2 :

 12 Fig. 1: The 3D-PCC. (A to C) Schematic diagrams for polariton frequencies as a function of detuning for a cQED system with two cavity modes (i.e., N = 2, red dashed lines) and one matter excitation (blue dashed line). When the light-matter coupling strength, g, is smaller than the frequency spacing between the cavity modes (∆ω) (A), the vacuum Rabi splittings are well separated in frequency. When g becomes larger than ∆ω (the superstrong coupling regime), there are two scenarios for the polariton branches, depending on the spatial overlap of the mode profiles: (B) the 2N scenario; and (C) the N + 1 scenario. UP, upper polariton; LP, lower polariton; MP, middle polariton. (D) Schematic of a 2DEG embedded in a 3D woodpile cavity. (E) Projected band structure of the bare woodpile cavity ((D) without the 2DEG layer). The grey-shaded regions and black lines denote modes extended in the whole crystal and the guided modes within the band gap, respectively. (F) Transmittance spectrum of the bare 3D-PCC from numerical simulation. The shaded region represents the photonic band gap expected from (E). The blue arrows mark cavity modes 1 to 4. The inset is a magnified view of the spectrum. (G) A photo of a fabricated architecture containing the designed woodpile structure.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Polarization-dependent mixing of photonic modes. (A to F) Magnetic-fielddependent transmittance (T ) spectra obtained from (A and B) experiments and calculations using (C and D) the full Hamiltonian and (E and F) the decoupled Hamiltonian, respectively. The top (bottom) panels are for TM (TE) modes. The white dots denote the peak frequencies extracted from experimental data using longer time-domain traces (24). The white dashed line shows the bare CR frequency. The LP-UP splitting and S-shaped MP are marked. The yellow arrows in (D) and (F) indicate the UP that shows a significant deviation between the spectra calculated using the full and decoupled Hamiltonian.

Contents 1 and Methods 3 1. 1

 131 Materials Preparation of the multiple quantum well sample . . . .

  The near-infrared (775 nm) output beam of a Ti:sapphire regenerative amplifier (1 kHz, 200 fs, Clark-MXR, Inc.,CPA-2001) was split into pump and probe beams. The pump beam was used to generate linearly polarized THz radiation in a ZnTe crystal through optical rectification. The generated THz beam was focused onto the sample that was mounted on the sample holder of a 10-T superconducting magnet cryostat (Oxford Instruments, Inc.). The direction of the static magnetic field was normal to the sample surface. The time-domain waveform of the transmitted THz radiation, E sample (t), was measured with another ZnTe crystal via electro-optic sampling with controlled time delays. A reference signal, E reference (t), was measured by repeating the measurements in the absence of a sample. E sample (t) and E reference (t) were Fourier-transformed into complex-valued frequencydomain spectra, Ẽsample (ω) and Ẽreference (ω), respectively. Transmittance, T , is equal to| Ẽsample (ω)| 2 /| Ẽreference (ω)| 2 .S2 Supplementary TextS2.1 Transmittance spectrum for the bare 3D-PCCWe performed transmission measurements on a bare 3D-PCC, i.e., a bare GaAs substrate without the 2DEG layer, to examine the quality of the fabricated 3D-PCC, see Fig.S1. The GaAs substrate was polished down to ∼85 µm due to the inaccuracy of the polishing process. Although the bare GaAs substrate has a thickness that was slightly different from the actual GaAs sample with the 2DEG layer, the thickness discrepancy only caused shifts of photonic mode frequencies. The frequencies of modes 1 and 2 in the experimental data are consistent with the simulation with an 85-µm-thick GaAs layer. Modes 3 and 4, which have extremely nar-row linewidths, were not observed in the experimental data due to the limited resolution of THz spectroscopy. The experimental data confirmed that the transmittance spectra of the bare 3D-PCC for TE and TM modes are degenerate.

ξ

  (z) is a normalized wave function that describes the QW confinement in the z direction. c † n,k and c n,k are the creation and annihilation operators of an electron in the state (n, k).

  iGxkl 2 c e -iG•ρ , which promotes an electron with momentum k in the highest-occupied Landau level (LL) n = ν -1 to the lowest-unoccupied LL n = ν and momentum k -G y . ν = 2πn e l 2 c denotes the filling factor of the 2DEG and N is the LL degeneracy. The in-plane position vector ρ is restricted to a woodpile unit cell. The CR excitation operators b(ρ) and b † (ρ) satisfy bosonic commutation relations when the number of CR excitations remains small compared to N , i.e., ⟨[ b(ρ), b † (ρ ′ )]⟩ = δ(ρρ ′ ), where ⟨• • • ⟩ denotes the expectation value in the electronic ground state (in which the lowest ν LLs are fully occupied). Note that the LL degeneracy entering the definition of the CR excitation operators is here written as N = a 2 /(2πl 2 c ), using the lattice constant a instead of the length L of the 2DEG. This is indeed the only physical choice as the light-matter coupling strength should not depend explicitly on any length scale other than l c in the plane. With those definitions, one recovers the light-matter coupling term H int discussed in the main text. The A 2 term is derived in a similar fashion.

)

  While we keep the same parameters (n e , ω p , ...) as in the model describing the woodpile structure, the new parameter ϵ ∈ [0, 1] allows us to artificially tune the spatial overlap between the two cavity modes. As ϵ is increased from 0 (no overlap) to 1 (perfect overlap), the splitting between the UP of the first mode p = 1 and the LP of the second one p = 2 becomes narrower and vanishes at ϵ = 1. The B-dependent transmission spectrum obtained with perfect overlap (ϵ = 1) is shown in S2A as an example. In this case, the polariton dispersion features a S-shaped MP in between the LP and the UP. The spectra found for ϵ ̸ = 1 are very similar to the ones reported in the main text with the full model describing the 3D-PCC. S2B shows the scaling of the intermode correlations ⟨a † 1 a 2 ⟩ computed from the toy model with the off-diagonal FOMη 12 = (dρ/a 2 ) g 1 (ρ)g 2 (ρ) ω c (ω 1 + ω 2 )/2introduced in the main text. The colored dots are obtained by tuning the mode frequencies (for ϵ = 1), while the black dashed line is obtained by tuning the overlap parameter ϵ. We observe that close to resonance, ω c ≈ ω p , intermode correlations are solely governed by η 12 , which is thus a relevant FOM for multimode ultrastrong photon-photon coupling.

  with the electric field operatorÊ(ρ, z) = i p,j ℏω p 2ε 0 a 3 E p,j (ρ, z)e j âp,jâ † p,j ,(S7)and e j the unit vector in the direction j = x, y. The contribution of the cavity mode p to the variance of the electric energy density [i.e., I(z) = p I p (z)] thus readsI p (z) = ℏω p 4a dρ a 2 ε(ρ, z) E 2 p,x (ρ, z) + E 2 p,y (ρ, z) . (S8) Note that the integral runs over a unit cell of the woodpile structure. The standard deviation of the in-plane spatial profile of the electric field (at the vertical location of the 2DEG) in the vacuum state reads E(ρ) = ⟨0| Ê2 (ρ, z 2DEG ) |0⟩. Similarly as before, one can write the contribution of the cavity mode p to E(ρ) [i.e., E(ρ) = p E p (ρ)] as

  standard deviation of the electric energy density along zσ p,j = dρ a 2 dz a (zzp,j ) 2 ε(ρ, z)E 2 p,j (ρ, z), ρ, z)E 2 p,j (ρ, z). (S12)Note that since photonic modes are localized in the vicinity of the defect layer in the z direction and exhibit discrete translational symmetry in the plane, we compute the spatial profile E p,j (ρ, z) in a woodpile unit cell of surface a 2 = (333 µm)2 . The integral along z runs over the whole computational cell, which includes air layers on each side of the woodpile structure in the z direction. The standard deviations of the 3D-PCC are σ 1,x ≈ σ 2,x ≈ σ 1,y ≈ σ 2,y = 0.6a

Fig. S1 :

 S1 Fig.S1: Transmittance spectrum for the bare 3D-PCC with a 85-µm-thick GaAs layer. We performed transmission measurements on a bare 3D-PCC, i.e., the defect layer is a bare GaAs substrate without the 2DEG layer. The red, black, and green lines denote the spectra obtained from numerical simulation and experiments with y-(TE) and x-(TM) polarized incident THz radiation, respectively. The blue arrows mark cavity modes 1 to 4.

Fig. S2 :Fig. S5 :

 S2S5 Fig. S2: Transmittance spectra and ground-state correlation of a toy model with two cavity modes and continuously tunable ξ p,p ′ . (A) Magnetic-field-dependent transmittance spectrum with perfect overlap (ϵ = 1), which shows an S-shaped MP. (B) Scaling of the intermode correlations with the off-diagonal coupling strength η 12 . The green and blue dots are obtained by tuning the frequency ω 1 = 2π × [0.339, 90] THz (while keeping ω 2 = 2π × 0.384 THz) and ω 2 = 2π × [0.384, 90] THz (while keeping ω 1 = 2π × 0.339 THz), respectively. The overlap parameter is ϵ = 1. The black dashed line is obtained by tuning ϵ between 0 and 1, with ω 1 = 2π × 0.339 THz and ω 2 = 2π × 0.384 THz. The magnetic field is set to B = 0.81 T, for which ω c = ω 1 . The overlap of all traces confirms that the intermode correlations at the ground state are governed by the FOM η 12 .

Fig. S6 :Fig. S7 :

 S6S7 Fig. S6: Simulations with different lengths of time window. (A) Transmittance spectra of a bare 3D-PCC with a series of time windows obtained from simulations. The peaks corresponding to the third and fourth cavity modes become pronounced when the time window is longer than 1000 ps. (B) Transmittance spectra in (A) that are plotted in logarithmic scale.

  Ga 0.76 As barriers. Silicon dopants were placed 80 nm away from the GaAs MQW. An electron density per well of 3.08 × 10 11 cm -2 and a mobility of 2.36 × 10 7 cm 2 /Vs were extracted from Hall transport measurements at 300 mK in the dark. The total electron density of the MQW was 3.08 × 10 12 cm -2 .

	S1 Materials and Methods
	S1.1 Preparation of the multiple quantum well sample
	. . . . . . . . . . . . A wafer containing multiple GaAs quantum wells (MQW) was grown by the Purdue molecular 3
	1.2 Fabrication of the THz woodpile cavity . . . . . . . . . . . . . . . . . . . . . beam epitaxy system. This structure had ten 30-nm-thick GaAs MQW separated by 160-nm 3
	1.3 THz time-domain magnetospectroscopy measurements . . . . . . . . . . . . . Al 0.24	4
	2 Supplementary Text	4
	3 Fig. S1 to S8	15

2.1 Transmittance spectrum for the bare 3D-PCC . . . . . . . . . . . . . . . . . . 4 2.2 Microscopic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 Ground-state correlations from different terms in Hamiltonian . . . . . . . . . 9 2.4 Vacuum fluctuations and localization of the cavity modes . . . . . . . . . . . . 10 2.5 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.6 Thickness of the defect layer in numerical calculations . . . . . . . . . . . . . 12 2.7 Simulations with different ranges of time delays . . . . . . . . . . . . . . . . . 12 2.8 Extraction of peak frequencies . . . . . . . . . . . . . . . . . . . . . . . . . . 13

Acknowledgments

This work was done in part using resources of the Research Support Shop and Shared Equipment Authority at Rice University, as well as computational resources of the Centre de calcul de l'université de Strasbourg (CCUS). We thank Motoaki Bamba, Kaden Hazzard, Tal Schwartz, Thibault Chervy, and Cyriaque Genet for useful discussions.

Funding: J.K. acknowledges support from the U.S. Army Research Office (through Award No. W911NF2110157), the Gordon and Betty Moore Foundation (through Grant No. 11520), the W. M. Keck Foundation (through Award No. 995764), and the Robert A. Welch Foundation (through Grant No. C-1509).

Authors' contributions: F.T., D.H., and J.K. conceptualized the project. F.T. and A.M. designed and fabricated the cavity device. F.T., A.B., and H.X. performed the measurements. F.T. analyzed the experimental data. D.H. derived the microscopic model and performed MEEP simulations and calculations. S.L., G.C.G., and M.J.M. grew the 2DEG sample by the molecular beam epitaxy system. F.T., S.S., and A.A. conducted COMSOL simulations. A.B., D.H., and J.K. supervised the project. F.T., D.H., and J.K. wrote the manuscript, with inputs from all authors.

Competing interests: The authors have no competing interests.

Data and materials availability: All data and codes are deposited at Dryad.

Supplementary Materials

Materials and Methods

Supplementary Text


R Thanikonda
email: ravi.thanikonda@student.unisi.it

J Ruiz-Garcia

D González-Ovejero

G Toso

E Martini

S Maci

Low Profile Wideband Scanning Antenna in Ka-band Based on Fully Metallic Metasurfaces

In this paper, a high gain antenna based on low profile metasurfaces is designed and numerically analyzed. The operation principle is the same of series fed Continuous Transverse Stub (CTS) structures, but in comparison to prior designs, radiation efficiency and compactness are improved by feeding the radiating aperture with a planar lens working in reflection. The lens can be realized using a pin-type metasurface, resulting in a fully metallic construction. The radiating aperture is realized by etching slots in the upper wall of a parallel plate waveguide loaded with corrugations; this solution provides an almost frequency-and scan-independent active impedance, resulting in wideband performance for all the pointing directions. Furthermore, because of the axial symmetry of the feeding structure, the antenna has the capability of generating multiple beams from the same aperture while maintaining excellent port decoupling.

I. INTRODUCTION

Satellite mobile communications have grown at an exponential rate in recent years as a result of the introduction of numerous new mega-constellations. In the ground segment, this calls for low-cost high-gain antennas with wide beam-scanning capabilities, as well as a low profile and low power consumption. From a technological standpoint, phased array antennas might be an option; nevertheless, they are currently too expensive and power-demanding for the consumer market. Traditional mechanically scanned reflectors, on the other hand, are too large and heavy for mobile platforms.

Recently, low profile solutions based on unconventional mechanical scanning mechanisms have been proposed [START_REF] Tekkouk | Wideband and large coverage continuous beam steering antenna in the 60-GHz band[END_REF]. In this solution, a radiating aperture consisting of a slotted parallel plate waveguide is excited by a corporate feeding network sitting on an underlying layer. Through a relative or common rotation of the two layers the beam can be scanned continuously in the upper hemisphere with no need for phase shifters.

A novel solution based on this concept is presented in this paper. The main novelty is relevant to the replacement of the corporate feeding network with a metasurface-based lens working in reflection, according to the principle presented in [2]. This solution allows for an extreme compactness of the overall antenna structure; furthermore, thanks to the axial symmetry of the lens, it provides multibeam capability. The complete structure consists of two parts: a Reflecting Luneburg Lens (RLL) realized in the bottom layer serves as a feeding structure, while a periodic array of slots is placed in the upper layer as a radiating structure. The two parts are coupled by a corner reflector.

Radiating Slot

Corrugations

A. Feeding Structure

The radiating part on the top is fed by a lens placed in the bottom layer, designed in circular shape and realized in a parallel plate waveguide (PPW). The lens is designed according to the operation principle introduced in [START_REF] Tekkouk | Wideband and large coverage continuous beam steering antenna in the 60-GHz band[END_REF]: the azimuthally symmetric modulation of the effective refractive index is designed to collimate in the upper layer, after reflection at the lens boundary, the rays launched by a point source located at the periphery (see Fig. 2). A more practical and compact structure is obtained by using as primary feed a printed horn positioned at the internal quasifocus point. The modulation of the refractive index profile can be obtained by loading one wall of the parallel plate waveguide with a spatially variable metasurface. This latter can be realized either by printing electrically small patches over a dielectric slab or using metallic pins. The latter solution has the advantage to lead to a fully metallic structure, thus, reducing losses.

B. Radiating Structure

The radiating structure sitting in the upper layer is also realized in PPW, and can be developed by hollow waveguide technology. This PPW is loaded with corrugations to support a slow wave. It is excited with a planar wavefront by the underlying feeding structure, through a corner reflector. The relative rotational movement of the radiating and feeding structures, described by the angle ϕ, creates a variable phase gradient along the radiating slots, which is responsible for the beam scanning. Fig. 3 shows the geometry of the unit cell of the radiating structure, which contains a reflection cancelling notch to allow for broadside radiation when the period is equal to the guided wavelength. By adjusting the parameters h notch & u we can easily control the reflection coefficient of the unit-under all scanning conditions and in a broadband. The antenna structure has been designed with overall radius of 7λ and radiating aperture radius 6.1909λ at 29GHz. At the design frequency, the reflection coefficient is lower than -15dB and the coupling coefficient between facing ports is -15dB. Fig. 4 presents the simulated radiation patters under scanning conditions in the range between -60 • to +60 • , corresponding to rotation angles from ϕ = -40 • to ϕ = +40 • . within the scanning range the overall sidelobes are lessthan 17dB. And the total efficiency is around 67% at broadside.

Fig. 1 :

 1 Fig. 1: 3D view of the complete design

Fig. 2 :

 2 Fig. 2: Reflecting Luneburg Lens (RLL).

Fig. 3 :

 3 Fig. 3: Reflection cancelling unit cell with corrugated PPW. Unit-cell parameters are h cppw = 2.5mm, w slot = 4.6mm, h slot = 2mm, h notch = 1mm

  performance in the E-plane at f = 29GHz. The step size of rotation angle is ∆ϕ = 10 • IV. CONCLUSION A compact, low profile high gain scanning antenna has been designed and numerically verified. Scanning is based on the series CTS principle, here combined with an innovative and effective solution for the beamforming network, based on a metasurface-based reflecting lens. ACKNOWLEDGEMENT This work was supported by Italian and European Space Agencies through contract 000133153/20/NL/NR.


Hamed Kabiri Kenari
email: hamedkabiri@ieee.org

Arash Ahadi
email: arash.ahadi@nimasystem.com

Independent Reaserchers

Multilayer Heterogeneous Graph Attention Network

Keywords: Heterogeneous Network, Neural Network, Graph Analysis

Graph neural network, a powerful graph representation technique based on deep learning, has shown superior performance and attracted considerable research interest. It has been considered in graph neural networks for a heterogeneous graph, which contains different types of nodes and links. However, these graph neural networks are based on choosing meta-paths. When the heterogeneous graphs have a lot of node types and edge types, choosing meta-paths is difficult. One of the most exciting advancements in deep learning is the attention mechanism, whose great potential has been well demonstrated in various areas. In this paper, we look at heterogeneous networks as multilayer networks, which means each edge type shows one layer. We propose a novel heterogeneous graph neural network based on hierarchical attention, including intra-layer and inter-layer aggregations. Specifically, intra-layer aggregation aims to learn the importance between a node and its layer-based neighbours, while inter-layer aggregation is able to learn the importance of different layers. With the learned importance from both intra-layer and inter-layer attention, the importance of node and layer can be fully considered. Then, the proposed model can generate node embedding by aggregating features from layer-based neighbours in a hierarchical manner. Extensive experimental results on two real-world heterogeneous graphs show the superior performance of our proposed model over the state-of-the-art.

INTRODUCTION

Many real-world datasets are naturally represented in a graph data structure, where objects and the relationships among them are embodied by nodes and edges, respectively. Examples include social networks [START_REF] Hamilton | Inductive Representation Learning on Large Graphs[END_REF][START_REF] Wang | Structural Deep Network Embedding[END_REF], physical systems [START_REF] Battaglia | Interaction Networks for Learning about Objects, Relations and Physics[END_REF][START_REF] Fout | Protein Interface Prediction using Graph Convolutional Networks[END_REF], traffic networks [START_REF] Li | Diffusion Convolutional Recurrent Neural Network: Data-Driven Traffic Forecasting[END_REF][START_REF] Zhang | GaAN: Gated Attention Networks for Learning on Large and Spatiotemporal Graphs[END_REF], citation networks [START_REF] Atwood | Diffusion-Convolutional Neural Networks[END_REF][START_REF] Hamilton | Inductive Representation Learning on Large Graphs[END_REF][START_REF] Kipf | Semi-Supervised Classification with Graph Convolutional Networks[END_REF], recommender systems [START_REF] Van Den | Graph Convolutional Matrix Completion[END_REF][START_REF] Zhang | STAR-GCN: Stacked and Reconstructed Graph Convolutional Networks for Recommender Systems[END_REF], knowledge graphs [START_REF] Bordes | Translating Embeddings for Modeling Multi-relational Data[END_REF][START_REF] Sun | RotatE: Knowledge Graph Embedding by Relational Rotation in Complex Space[END_REF], and so on. The unique non-Euclidean nature of graphs renders them difficult to be modelled by traditional machine learning models. There is no order or size limit for each node's neighbourhood set. However, most statistical models assume an ordered and fixed-size input in the Euclidean space. Therefore, nodes could be represented by meaningful low-dimensional vectors in Euclidean space and then be taken as the input for other machine learning models.

Different graph embedding techniques have been proposed for the graph structure. LINE [START_REF] Tang | LINE: Largescale Information Network Embedding[END_REF] generates node embeddings by exploiting the first-order and second-order proximity between nodes. Random-walkbased methods including DeepWalk [START_REF] Perozzi | DeepWalk: Online Learning of Social Representations[END_REF], node2vec [START_REF] Grover | Node2Vec: Scalable Feature Learning for Networks[END_REF], and TADW [START_REF] Yang | Network Representation Learning with Rich Text Information[END_REF] feed node sequences generated by random walks to a skip-gram model [START_REF] Mikolov | Efficient Estimation of Word Representations in Vector Space[END_REF] to learn node embeddings. With the rapid development of deep learning, graph neural networks (GNNs) have been proposed, which learn the graph representations using specially designed neural layers. Spectral-based GNNs, including ChebNet [START_REF] Defferrard | Convolutional Neural Networks on Graphs with Fast Localized Spectral Filtering[END_REF] and GCN [START_REF] Kipf | Semi-Supervised Classification with Graph Convolutional Networks[END_REF], perform graph convolution operations in the Fourier domain of an entire graph. Recent spatial-based GNNs, including GraphSAGE [START_REF] Hamilton | Inductive Representation Learning on Large Graphs[END_REF], GAT [START_REF] Velickovic | Graph Attention Networks[END_REF], and many other variants [START_REF] Li | Gated Graph Sequence Neural Networks[END_REF][START_REF] Zhang | GaAN: Gated Attention Networks for Learning on Large and Spatiotemporal Graphs[END_REF][START_REF] Zhang | STAR-GCN: Stacked and Reconstructed Graph Convolutional Networks for Recommender Systems[END_REF], address the scalability and generalization ability of the spectral-based models by performing graph convolution operations directly in the graph domain. An increasing number of researchers have paid attention to this promising area.

Although GNNs have achieved state-of-the-art results in many tasks, most GNN-based models assume that the input is a homogeneous graph with only one node type and one edge type. Most real-world graphs consist of various types of nodes and edges associated with attributes in different feature spaces. For example, a co-authorship network contains at least two types of nodes: authors and papers. Author attributes may include affiliations, citations, and research fields. Paper attributes may include keywords, venue, year, etc. We refer to graphs of this kind as heterogeneous information networks (HINs) or heterogeneous graphs. The heterogeneity in graph structure and node content makes it challenging for GNNs to encode their rich and diverse information into a low-dimensional vector space.

Most existing heterogeneous graph embedding methods are based on the idea of metapaths. A metapath is an ordered sequence of node types and edge types defined on the network schema, which describes a composite relation between the node types involved. For example, in a scholar network with authors, papers, and venues, Author-Paper-Author (APA) and Author-Paper-Venue-Paper-Author (APVPA) are metapaths describing two different relations among authors. The APA metapath associates two co-authors, while the APVPA metapath associates two authors who published papers in the same venue. Therefore, we can view a metapath as a high-order proximity between two nodes. Because traditional GNNs treat all nodes equally, they cannot model the complex structural and semantic information in heterogeneous graphs.

Although these metapath-based embedding methods outperform traditional network embedding methods on various tasks, such as node classification and clustering, they still suffer from at least one of the following limitations. (1) The model does not leverage node content features, so it rarely performs well on heterogeneous graphs with rich node content features (e.g., metapath2vec [START_REF] Dong | Metapath2Vec: Scalable Representation Learning for Heterogeneous Networks[END_REF], ESim [START_REF] Shang | Meta-Path Guided Embedding for Similarity Search in Large-Scale Heterogeneous Information Networks[END_REF], HIN2vec [START_REF] Fu | HIN2Vec: Explore Meta-paths in Heterogeneous Information Networks for Representation Learning[END_REF], and HERec [START_REF] Shi | Heterogeneous Information Network Embedding for Recommendation[END_REF]). [START_REF] Wang | Structural Deep Network Embedding[END_REF] The model requires a manual metapath selection process when number of node types and edges types are high, choosing methapaths is difficult (e.g., metapath2vec [START_REF] Dong | Metapath2Vec: Scalable Representation Learning for Heterogeneous Networks[END_REF], HERec [START_REF] Shi | Heterogeneous Information Network Embedding for Recommendation[END_REF] and HAN [START_REF] Wang | Heterogeneous Graph Attention Network[END_REF])).

To address these limitations, we propose a novel Metapath Aggregated Graph Neural Network (MLHAN) for heterogeneous graph embedding. MLHAN addresses all the issues described above by applying node content transformation, intra-layer aggregation, and inter-layer aggregation to generate node embeddings. Specifically, MLHAN first applies type-specific linear transformations to project heterogeneous node attributes, with possibly unequal dimensions for different node types, to the same latent vector space. Next, MLHAN applies intra-layer aggregation with the attention mechanism [START_REF] Velickovic | Graph Attention Networks[END_REF] for every layer. During this intra-layer aggregation, each target node extracts and combines information from the node with its layerbased neighbors. In this way, MLHAN captures the structural and semantic information of heterogeneous graphs from both neighbor nodes. Following intra-layer aggregation, MLHAN further conducts inter-layer aggregation using the attention mechanism to fuse latent vectors obtained from multiple layers into final node embeddings. By integrating multiple layers, our model can learn the comprehensive semantics ingrained in the heterogeneous graph. In summary, this work makes several major contributions:

(1) We propose a novel multilayer aggregated graph neural network for heterogeneous graph embedding.

(2) We conduct extensive experiments on the IMDb and the DBLP datasets for node classification and node clustering to evaluate the performance of our proposed model. Experiments on all of these datasets and tasks show that the node embeddings learned by MLHAN are consistently better than those generated by other state-of-the-art baselines.

Graph Neural Network

Graph neural networks (GNNs) which aim to extend the deep neural network to deal with arbitrary graphstructured data are introduced in [START_REF] Gori | A new model for learning in graph domains[END_REF][START_REF] Scarselli | The graph neural network model[END_REF]. [START_REF] Li | Gated Graph Sequence Neural Networks[END_REF] proposes a propagation model incorporating gated recurrent units to propagate information across all nodes. Recently, there has been a surge of generalizing convolutional operations on graph-structured data. The graph convolutional neural work generally falls into two categories: spectral domain and non-spectral domain. On one hand, spectral approaches work with a spectral representation of the graphs. In [START_REF] Bruna | Spectral networks and locally connected networks on graphs[END_REF] extends convolution to general graphs by finding the corresponding Fourier basis. In [START_REF] Defferrard | Convolutional Neural Networks on Graphs with Fast Localized Spectral Filtering[END_REF] utilizes K-order Chebyshev polynomials to approximate smooth filters in the spectral domain. [START_REF] Kipf | Semi-Supervised Classification with Graph Convolutional Networks[END_REF] proposes a spectral approach named Graph Convolutional Network, which designs a graph convolutional network via a localized first-order approximation of spectral graph convolutions. On the other hand, we also have non-spectral approaches, which define convolutions directly on the graph, operating on groups of spatially close neighbors. [START_REF] Hamilton | Inductive Representation Learning on Large Graphs[END_REF] introduces GraphSAGE, which performs a neural network-based aggregator over a fixed-size node neighbor. It can learn a function that generates embeddings by aggregating features from a node's local neighborhood.

Attention mechanisms, e.g., self-attention [START_REF] Vaswani | Attention is All you Need[END_REF] and soft-attention [START_REF] Bahdanau | Neural machine translation by jointly learning to align and translate[END_REF], have become one of the most influential mechanisms in deep learning. Some previous works introduce the attention mechanism for graph based applications, e.g., the recommendation [START_REF] Han | Aspect-Level Deep Collaborative Filtering via Heterogeneous Information Networks[END_REF][START_REF] Hu | Leveraging Meta-path based Context for Top-N Recommendation with A Neural Co-Attention Model[END_REF]. Inspired by the attention mechanism, the Graph Attention Network [START_REF] Velickovic | Graph Attention Networks[END_REF] is proposed to learn the importance between nodes and their neighbors and fuse the neighbors to perform node classification. However, the above graph neural network cannot deal with various types of nodes and edges and can only be applied to homogeneous graphs.

Network Embedding

Network embedding, i.e., network representation learning (NRL), is proposed to embed a network into a low-dimensional space while preserving the network structure and property to apply the learned embeddings to the downstream network tasks. For example, the random walk based methods [START_REF] Perozzi | DeepWalk: Online Learning of Social Representations[END_REF][START_REF] Grover | Node2Vec: Scalable Feature Learning for Networks[END_REF], the deep neural network based methods [START_REF] Wang | Structural Deep Network Embedding[END_REF], the matrix factorization based methods [START_REF] Ou | Asymmetric transitivity preserving graph embedding[END_REF][START_REF] Xiaowang | Community Preserving Network Embedding[END_REF], and others, e.g., LINE [START_REF] Tang | LINE: Largescale Information Network Embedding[END_REF]. However, all these algorithms are proposed for the homogeneous graphs. Some elaborate reviews can be found in [START_REF] Cui | A survey on network embedding[END_REF][START_REF] Goyal | Graph embedding techniques, applications, and performance: A survey[END_REF].

Heterogeneous graph embedding mainly focuses on preserving the meta-path based structural information. ESim [START_REF] Shang | Meta-Path Guided Embedding for Similarity Search in Large-Scale Heterogeneous Information Networks[END_REF] accepts userdefined meta-paths as guidance to learn vertex vectors in a userpreferred embedding space for similarity search. Even through ESim can utilize multiple meta-paths, it cannot learn the importance of meta-paths. To achieve the best performance, ESim must conduct a grid search to find the optimal meta-path weights. It is hard to find the optimal combination for a specific task. Metapath2vec [START_REF] Dong | Metapath2Vec: Scalable Representation Learning for Heterogeneous Networks[END_REF] designs a meta-path-based random walk and utilizes skip-gram to perform heterogeneous graph embedding. However, metapath2vec can only utilize one meta-path and may ignore some helpful information. Like metapath2vec, HERec [START_REF] Shi | Heterogeneous Information Network Embedding for Recommendation[END_REF] proposes a type constraint strategy to filter the node sequence and capture the complex semantics reflected in a heterogeneous graph. HIN2Vec [START_REF] Fu | HIN2Vec: Explore Meta-paths in Heterogeneous Information Networks for Representation Learning[END_REF] carries out multiple prediction training tasks, which simultaneously learn the latent vectors of nodes and meta-paths. [START_REF] Chen | PME: Projected Metric Embedding on Heterogeneous Networks for Link Prediction[END_REF] proposes a projected metric embedding model, PME, which can preserve node proximities via Euclidian Distance. PME projects different types of nodes into the same relation space and conducts heterogeneous link prediction. To study the problem of comprehensively describing heterogeneous graphs, [START_REF] Shi | Easing Embedding Learning by Comprehensive Transcription of Heterogeneous Information Networks[END_REF] proposes HEER, which can embed heterogeneous graphs via edge representations. [START_REF] Fan | Gotcha-sly malware!: Scorpion a metagraph2vec based malware detection system[END_REF] proposes an embedding model metagraph2vec, where the structures and semantics are maximally preserved for malware detection.

In [START_REF] Sun | Joint embedding of meta-path and meta-graph for heterogeneous information networks[END_REF] proposes meta-graph-based network embedding models, which simultaneously consider the hidden relations of all meta information of a meta-graph.

PRELIMINARY

A heterogeneous graph is a special kind of information network, which contains either multiple types of objects or multiple types of links. Definition 3.1. Heterogeneous Graph [START_REF] Sun | Mining heterogeneous information networks: a structural analysis approach[END_REF]. A heterogeneous graph, denoted as 𝐺 = (𝑉, 𝐸), consists of an object set 𝑉 and a link set 𝐸. A heterogeneous graph is also associated with a node type mapping function 𝜙: 𝑉 → 𝐴 and a link type mapping function 𝜑: 𝐸 → 𝑅. 𝐴

THE PROPOSED MODEL

This section proposes a novel semi-supervised graph neural network for heterogeneous graphs. Our model follows a hierarchical attention structure: intra-layer aggregation → inter-layer aggregation.

First, we propose intra-layer attention to learn the weight of layer based neighbors and aggregate them to get the layer-specific node embedding. After that, MLHAN can tell the difference of layers via inter-layer attention and get the optimal weighted combination of the inter-layer node embedding for the specific task.

Intra-layer Aggregation

Before aggregating the information from intra-layer neighbors for each node, we should notice that the layer-based neighbors of each node play a different role and show different importance in learning node embedding for the specific task. Here, we introduce intra-layer attention to learn the importance of layerbased neighbors for each node in a heterogeneous graph and aggregate the representation of these meaningful neighbors to form a node embedding.

According to our definition of multilayer heterogeneous graph, each node in heterogeneous graph can has same node type neighbors in unipartite layers or different node type neighbors in bipartite layers.

Due to the heterogeneity of nodes, different types of nodes have different feature spaces. Therefore, for each type of node (e.g., node with type 𝐴), we design the type-specific transformation matrix 𝑀 𝐴 to project the features of different types of nodes into the same feature space. Unlike [START_REF] Hamilton | Embedding logical queries on knowledge graphs[END_REF], the type-specific transformation matrix is based on node type rather than edge type. The projection process can be shown as follows:

ℎ 𝑖 / = 𝑀 𝐴 . ℎ 𝑖 (1)
where ℎ 𝑖 and ℎ 𝑖 / are the original and projected features of the node 𝑖, respectively. By type-specific projection operation, the inter-layer attention can handle arbitrary types of nodes.

After that, we leverage self-attention [START_REF] Vaswani | Attention is All you Need[END_REF] to learn the weight among various kinds of nodes. Given a node pair (𝑖, 𝑗) which are connected via layer Φ, the intra-layer attention 𝑒 𝑖𝑗 Φ can learn the importance 𝑒 𝑖𝑗 Φ , which means how important the node 𝑗 will be for the node 𝑖. The importance of layer-based node pairs (𝑖, 𝑗) can be formulated as follows:

𝑒 𝑖𝑗 Φ = 𝑎𝑡𝑡 𝑖𝑛𝑡𝑟𝑎-𝑙𝑎𝑦𝑒𝑟 (ℎ 𝑖 / , ℎ 𝑗 / ; Φ) (2)
Here 𝑎𝑡𝑡 𝑛𝑜𝑑𝑒 denotes the deep neural network, which performs the intra-layer attention. The above Eq. ( 2)

shows that the weight of layer-based node pairs (𝑖, 𝑗) depends on their features in a given layer Φ.

where 𝜎 denotes the activation function, || denotes the concatenate operation, and 𝑎 Φ is the intra-layer attention vector for the layer Φ. As we can see from Eq. ( 3), the weight coefficient of (𝑖, 𝑗) depends on their features. Also, please note that the weight coefficient 𝛼 𝑖𝑗 Φ is asymmetric, which means they make different contributions to each other, not only because of the concatenated order in the numerator but also because they have different neighbors, so the normalized term (denominator) will be quite different.

Then, the layer-based embedding of a node 𝑖 can be aggregated by the neighbor's projected features with the corresponding coefficients as follows:

𝑧 𝑖 Φ = 𝜎 (∑ 𝛼 𝑖𝑗 Φ . ℎ 𝑗 / 𝑗𝜖𝑁 𝑖 Φ ) (4)
where 𝑧 𝑖 Φ is the learned embedding of a node 𝑖 in the layer Φ. Every node embedding is aggregated by its neighbors. Since the attention weight 𝛼 𝑖𝑗 Φ is generated for a single layer, it is layer-specific.

Since heterogeneous graphs present scale-free property, the variance of graph data is relatively high. To tackle the above challenge, we extend intra-layer attention to multihead attention to make the training process more stable. Specifically, we repeat the intra-layer attention for 𝐾 times and concatenate the learned embeddings as the layer-specific embedding:

𝑧 𝑖 Φ = || 𝑘=1 𝐾 𝜎 (∑ 𝛼 𝑖𝑗 Φ . ℎ 𝑗 / 𝑗𝜖𝑁 𝑖 Φ ) (5)
To sum up, given the projected feature vectors ℎ 𝑢 / ∈ ℝ 𝑑 / ∀𝑢 ∈ 𝒱 that has node type 𝐴 ∈ 𝒜 and the set of layers {Φ 1 , … , Φ 𝑃 }, after feeding node features into intra-layer attention, we can obtain layer-specific vector representations of the target node 𝑢 ∈ 𝒱 𝐴 , denoted as {ℎ 𝑣 Φ 1 , ℎ 𝑣 Φ 2 , … , ℎ 𝑣 Φ 𝑝 }.

Inter-layer Aggregation

To learn a more comprehensive node embedding, we need to fuse layers. To address the challenge of layer fusion in a heterogeneous graph, we propose a novel inter-layer attention to automatically learn the importance of different layers and fuse them for the specific task. Now, we have |𝒱 𝐴 | sets of latent vectors for a node type 𝐴: {ℎ 𝑣 Φ 1 , ℎ 𝑣 Φ 2 , … , ℎ 𝑣 Φ 𝑝 } for 𝑣 ∈ 𝒱 𝐴 , where 𝑃 is the number of layers for type 𝐴.

Taking 𝑃 groups of layer-specific node embeddings for type 𝐴 learned from intra-layer attention as input, the learned weights of each layer (𝛽 𝐴 Φ 1 , … , 𝛽 𝐴 Φ 𝑝 ) can be shown as follows:

(𝛽 𝐴 Φ 1 , … , 𝛽 𝐴 Φ 𝑝 ) = 𝑎𝑡𝑡 𝑖𝑛𝑡𝑒𝑟-𝑙𝑎𝑦𝑒𝑟 (ℎ 𝐴 Φ 1 , ℎ 𝐴 Φ 2 , … , ℎ 𝐴 Φ 𝑝 ) (6)
Here 𝑎𝑡𝑡 𝑖𝑛𝑡𝑒𝑟-𝑙𝑎𝑦𝑒𝑟 denotes the deep neural network which performs the inter-layer attention. It shows that inter-layer attention can capture various types of relations information behind a heterogeneous graph.

To learn the importance of each layer, we first transform layer-specific embedding through a nonlinear transformation (e.g., one-layer MLP). Then, we measure the importance of layer-specific embedding as the similarity of transformed embedding with a type attention vector 𝑞 𝐴 . Furthermore, we average the importance of all the layer-specific node embedding, which can be explained as the importance of each layer for all nodes. The importance of each layer, denoted as 𝑤 Φ 𝑃 , is shown as follows:

𝑤 Φ 𝑝 = 1 |𝒱 𝐴 | ∑ (𝑞 𝐴 ) 𝑇 𝑖𝜖𝒱 𝐴 . tanh (𝑊 𝐴 . ℎ 𝑖 Φ 𝑝 + 𝑏 𝐴 ) (7)
Where 𝒱 𝐴 is nodes set in type 𝐴, 𝑊 𝐴 is the type specific weight matrix, 𝑏 𝐴 is the type specific bias vector, 𝑞 𝐴 is the type specific attention vector. After obtaining the importance of each layer, we normalize them via softmax function. The weight of layer Φ 𝑠 , denoted as 𝛽 Φ 𝑠 , can be obtained by normalizing the above importance of all layers using softmax function,

𝛽 Φ 𝑠 = 𝑒𝑥𝑝(𝑤 Φ 𝑠 ) ∑ 𝑒𝑥𝑝(𝑤 Φ 𝑟 ) 𝑃 𝑟=1 (8)
The higher 𝛽 Φ 𝑠 , the more important layer Φ 𝑠 is. With the learned weights as coefficients, we can fuse these layer-specific embeddings to obtain the final embedding ℎ 𝑣 as follows:

ℎ 𝑣 = ∑ 𝛽 Φ 𝑖 . ℎ 𝑣 Φ 𝑖 𝑃 𝑖=1 (9)
The final embedding is aggregated by all layer-specific embedding. After applying components introduced in the previous sections, we obtain the final node representations, which can then be used in different downstream tasks. We can train MLHAN in semi-supervised learning. For semi-supervised learning, with the guide of a small fraction of labeled nodes, we can optimize the model weights by minimizing the cross entropy via backpropagation and gradient descent, and thereby learn meaningful node embeddings for heterogeneous graphs. The cross entropy loss for this semi-supervised learning is formulated as: [START_REF] Van Den | Graph Convolutional Matrix Completion[END_REF] where 𝑉 𝐿 is the set of nodes that have labels, 𝐶 is the number of classes, 𝑦 𝑣 is the one-hot label vector of node 𝑣, and ℎ 𝑣 is the predicted probability vector of node 𝑣.

ℒ = -∑ ∑ 𝑦 𝑣 [𝑐]. log ℎ 𝑣 [𝑐] 𝐶 𝑐=1 𝑣∈𝑉 𝐿

EXPERIMENTS

Datasets

The detailed descriptions of the heterogeneous graph used here are shown in Table 1.

• DBLP is a computer science bibliography website. We adopt a subset of DBLP extracted by [START_REF] Gao | Graph-based Consensus Maximization Among Multiple Supervised and Unsupervised Models[END_REF][START_REF] Ji | Graph Regularized Transductive Classification on Heterogeneous Information Networks[END_REF], containing 4057 authors, 14328 papers, 7723 terms, and 20 publication venues after data preprocessing. The authors are divided into four research areas (Database, Data Mining, Artificial Intelligence, and Information Retrieval). Each author is described by a bag-of-words representation of their paper keywords.

• IMDB is an online database about movies and television programs, including information such as cast, production crew, and plot summaries. We use a subset of IMDb scraped from online, containing 4278 movies, 2081 directors, and 5257 actors after data preprocessing. Movies are labeled as one of three classes (Action, Comedy, and Drama) based on their genre information. Each movie is also described by a bag-ofwords representation of its plot keywords.

Baselines

We compare with some state-of-art baselines, including the (heterogeneous) network embedding methods and graph neural network-based methods, to verify the effectiveness of the proposed MLHAN.

• DeepWalk [START_REF] Perozzi | DeepWalk: Online Learning of Social Representations[END_REF]: A random walk-based network embedding method which designs for homogeneous graphs. Here, we ignore the heterogeneity of nodes and perform DeepWalk on the whole heterogeneous graph. . node2vec [START_REF] Grover | Node2Vec: Scalable Feature Learning for Networks[END_REF] is a traditional homogeneous model serving as a generalized version of DeepWalk [START_REF] Perozzi | DeepWalk: Online Learning of Social Representations[END_REF]. We apply it to the heterogeneous graphs by ignoring the heterogeneity of the graph structure and dropping all node content features.

• ESim [START_REF] Shang | Meta-Path Guided Embedding for Similarity Search in Large-Scale Heterogeneous Information Networks[END_REF]: A heterogeneous graph embedding method that can capture semantic information from multiple meta-paths. Because it is difficult to search the weights of a set of meta-paths, we assign the weights learned from HAN to ESim.

• metapath2vec [START_REF] Dong | Metapath2Vec: Scalable Representation Learning for Heterogeneous Networks[END_REF]: A heterogeneous graph embedding method that performs meta-path-based random walk and utilizes skip-gram to embed the heterogeneous graphs. Here, we test all the metapaths for metapath2vec and report the best performance.

• HERec [START_REF] Shi | Heterogeneous Information Network Embedding for Recommendation[END_REF]: A heterogeneous graph embedding method that designs a type constraint strategy to filter the node sequence and utilizes skip-gram to embed the heterogeneous graphs. Here, we test all the meta-paths for HERec and report the best performance.

• GCN [START_REF] Kipf | Semi-Supervised Classification with Graph Convolutional Networks[END_REF]: It is a semi-supervised graph convolutional network designed for homogeneous graphs. Here, we test all the meta-paths for GCN and report the best performance.

• GAT [START_REF] Velickovic | Graph Attention Networks[END_REF]: It is a semi-supervised neural network that considers the attention mechanism on the homogeneous graphs. Here, we test all the meta-paths for GAT and report the best performance.

• HAN [START_REF] Wang | Heterogeneous Graph Attention Network[END_REF]: The proposed semi-supervised graph neural network simultaneously employs nodelevel and semantic-level attention.

Implementation Details

For the proposed MLHAN, we randomly initialize parameters and optimize the model with Adam [START_REF] Diederik | Adam: A method for stochastic optimization[END_REF]. For the proposed MLHAN, we set the learning rate to 0.005, the regularization parameter to 0.001, the dimension of the inter-layer attention vector 𝑞 𝐴 to 128 for each type, the number of attention heads 𝐾 to 8, and the dropout of attention to 0.6. Moreover, we use early stopping with patience of 100, i.e., we stop training if the validation loss does not decrease for 100 consecutive epochs. For GCN and GAT, we optimize their parameters using the validation set. For semi-supervised graph neural networks, including GCN, GAT, HAN, and MLHAN, we split precisely the same training, validation, and test sets to ensure fairness. For random walk-based methods, including DeepWalk, ESim, metapath2vec, and HERec, we set the window size to 5, walk length to 100, walks per node to 40, and the number of negative samples to 5. We set the embedding dimension to 64 for all the above algorithms for a fair comparison.

Classification

Here, we employ the KNN classifier with 𝑘 = 5 to perform node classification. Since the variance of graphstructured data can be quite high, we repeat the process ten times and report the averaged Macro-F1 and Micro-F1 in Table 2 Based on Table 2, we can see that MLHAN achieves the best performance. ESim, which can leverage multiple meta-paths, performs better than metapath2vec for the traditional heterogeneous graph embedding method. Generally, graph neural network-based methods that combine the structure and feature information, e.g., GCN and GAT, usually perform better. To go deep into these methods, compared to simply average over node neighbors, e.g., GCN, HAN, and MLHAN. GAT, HAN, and MLHAN can weigh the information properly and improve the performance of the learned embedding. Compared to GAT, the proposed MLHAN, designed for heterogeneous graphs, shows its superiority. The above analysis shows that the proposed MLHAN achieves the best performance on all datasets. The results demonstrate that capturing the importance of nodes and layers in heterogeneous graph analysis is quite important.

Clustering

We also conduct the clustering task to evaluate the embeddings learned from the above algorithms. Once the proposed MLHAN is trained, we can get all the node embedding via feed-forward. Here, we utilize the KMeans to perform node clustering, and the number of clusters 𝐾 is set to the number of classes. We use the same ground truth as in node classification. Moreover, we adopt NMI and ARI to assess the quality of the clustering results. Since the performance of KMeans is affected by initial centroids, we repeat the process ten times and report the average results in Table 3. As can be seen in Table 3, we can find that MLHAN performs consistently much better than all baselines. Also, graph neural network-based algorithms usually achieve better performance. Based on the above analysis, the proposed MLHAN can give a comprehensive description of heterogeneous graph and achieve a significant improvements.

CONCLUSION

This paper tackles one fundamental problem in heterogeneous graph neural networks choosing meta-paths for specific problems. We propose a semi-supervised multilayer heterogeneous graph neural network based solely on the attention mechanism. The proposed MLHAN can capture complex structures behind heterogeneous graphs by looking at heterogeneous graphs as multilayer networks. The proposed model leverages intra-layer aggregation and inter-layer aggregation to learn the importance of nodes and layers.

Meanwhile, the proposed model utilizes the structural and feature information uniformly. Experimental results, including classification and clustering, demonstrate the effectiveness of MLHAN. Analyzing the proposed MLHAN has proven its good interpretability.

  and 𝑅 denote the sets of predefined object types and link types, where |𝐴| + |𝑅| > 2. We begin by defining a heterogeneous graph with 𝑀 different types of nodes as 𝐺 = (𝐺 𝑈 , 𝐺 𝐵 ) where 𝐺 𝑈 = {𝐿 𝑖 : 𝑖 ∈ {1,2, … , 𝑀}} is a set of 𝑀 unipartite graphs that each one interconnects one type of nodes. and 𝐺 𝐵 = {𝐿 𝑖𝑗 : 𝑖, 𝑗 ∈ {1,2, … , 𝑀}, 𝑖 ≠ 𝑗} is a set of bipartite graphs that each one interconnects two different types of nodes. Each layer 𝐿 𝑖 = (𝑉 𝑖 , 𝐸 𝑖 ) with 𝑉 𝑖 and 𝐸 𝑖 represent the nodes of type 𝑖 and inner edges between them, respectively. Similarly, we can have layers 𝐿 𝑖𝑗 = (𝑉 𝑖 , 𝑉 𝑖 , 𝐸 𝑖𝑗 ) where {𝐸 𝑖𝑗 ⊆ 𝑉 𝑖 × 𝑉 𝑗 : 𝑖, 𝑗 ∈ {1,2, … , 𝑀}, 𝑖 ≠ 𝑗} a bipartite graph that connects nodes of type 𝑖 and type 𝑗. We use the 𝑁 × 𝑁 adjacency matrix 𝐴 as adjacency matrix of 𝐺 that 𝑁 = ∑ 𝑀 𝑖=1 |𝑉 𝑖 | and the 𝑉 𝑚 × 𝑉 𝑙 adjacency matrix 𝐴 𝑚𝑙 that indicates edges of set 𝐸 𝑚𝑙 , which are related as

	𝐴 11 ⋯ 𝐴 1𝑀
	𝐴 = [ ⋮	⋱ ⋮	]
	𝐴 𝑀1 ⋯ 𝐴 𝑀𝑀

  Please note that 𝑒 𝑖𝑗Φ is asymmetric, i.e., the importance of node 𝑖 to node 𝑗 and the importance of node 𝑗 to node 𝑖 in layer Φ can be quite different. It shows that intra-layer attention can preserve the asymmetry, a critical property of heterogeneous graphs.

		𝑖𝑗 Φ via
	the softmax function:
	𝛼 𝑖𝑗 Φ = 𝑠𝑜𝑓𝑡𝑚𝑎𝑥 𝑖 (𝑒 𝑖𝑗 Φ ) =	𝑒𝑥𝑝(𝜎(𝑎 Φ 𝑇 . [ℎ 𝑖 / ||ℎ 𝑗 / ])) 𝑒𝑥𝑝(𝜎(𝑎 Φ 𝑇 . [ℎ 𝑖 / ||ℎ 𝑘 / ])) 𝑘𝜖𝑁 𝑖 ∑ Φ

Then, we inject the structural information into the model via masked attention, which means we only calculate the 𝑒 𝑖𝑗 Φ for nodes 𝑗𝜖 𝑁 𝑖 Φ , where 𝑁 𝑖 Φ denotes the layer based neighbors of node 𝑖 in layer Φ. After obtaining the layer-based node pairs' importance, we normalize them to get the weight coefficient 𝛼

Table 1 .

 1 Statistics of datasets.

.

  Table2. Experiment results (%) on the IMDb and DBLP datasets for the node classification task.

	Dataset	Metrics	Train%	DeepWalk	Node2vex	ESim	Metapath2vec	HErec	GCN	GAT	HAN	MLHAN
			20%	41.84	50.12	49.49	47.17	46.73	53.85	54.76	57.31	59.07
		Macro-F1	40% 60%	46.31 49.25	51.75 52.77	51.21 52.57	48.69 49.29	47.92 47.96	54.79 55.36	56.62 57.58	57.27 58.41	59.99 60.38
	IMDb		80% 20%	51.47 47.50	52.61 51.06	52.49 50.44	51.11 48.34	48.85 47.35	55.89 53.92	58.55 54.76	59.63 57.44	61.16 59.32
		Micro-F1	40% 60%	51.11 53.33	52.89 53.91	52.33 53.65	49.29 50.99	49.01 49.31	54.88 55.35	56.68 57.59	58.44 59.54	60.22 60.60
			80%	55.45	53.84	53.66	51.62	50.23	55.75	58.52	60.36	61.25
			20%	78.55	87.82	91.80	89.59	91.94	89.12	92.17	92.81	93.85
		Macro-F1	40% 60%	82.14 84.79	89.19 89.81	92.73 92.96	91.03 91.62	92.56 93.20	90.22 90.55	92.36 92.54	93.08 93.26	93.95 94.59
	DBLP		80% 20%	85.93 80.49	90.05 88.33	93.39 92.33	91.98 90.14	93.37 92.61	91.10 89.63	92.85 92.73	93.62 93.45	94.82 94.33
		Micro-F1	40% 60%	83.85 86.39	89.69 90.21	93.17 93.40	91.48 92.06	93.17 93.78	90.34 90.69	92.89 93.09	93.69 93.84	94.40 95.71
			80%	87.38	90.49	93.80	92.43	93.90	91.45	93.36	94.35	95.19

Table 3 .

 3 Experiment results (%) on the IMDb and DBLP datasets for the node clustering task.

	Dataset	Metrics	Node2vex	DeepWalk	ESim	Metapath2vec	HErec	GCN	GAT	HAN	MLHAN
		NMI	5.78	2.01	1.63	1.45	0.95	8.02	8.40	11.35	12.94
	IMDb
		ARI	6.58	2.71	1.57	0.78	0.67	8.25	9.43	11.67	14.10
		NMI	78.13	77.65	69.45	75.30	70.15	74.57	71.85	78.61	80.52
	DBLP
		ARI	82.49	92.47	73.34	79.23	73.57	78.62	77.16	84.07	85.26