Wave-vector and polarization dependent impedance model for a hexagonal periodic metasurface exemplified through finite-difference time-domain simulations

Yi S. Ding; Yang He

doi:10.1364/OE.25.020757

1. Introduction

Metasurfaces are thin screen structures composed of artificial unit cells of subwavelength size [1]. Planar periodic metallic surfaces (PMS) are usual cases of metasurfaces and they are usually termed as frequency selective surfaces for microwave frequencies [2]. Recently, these structures are playing more important roles in manipulating optical fields [3].

For infinite PMS, one can restrict theoretical analysis within the unit cell by use of the Floquet theorem. However, for finite large PMS in practice, the theoretical analysis (especially full-wave analysis) is computationally very expensive because of the large amount and fine details of unit cells [4]. One remedy is to use the high-frequency approximation, e.g. geometric and physics optics approximation [5]. These approximation can not consider multiple scattering from nearby scatterers. Another kind of approximation is to homogenize the PMS, just as what has been done for natural dielectric materials in macroscopic optics. Impedance sheet model (ISM), as described in [6, 7], is such a kind of approximation. The PMS is effectively homogenized as an infinitely thin sheet with an impedance defined as the ratio of the average tangential electric fields and average surface currents.

In practice, however, the surface impedance exhibits different values for different incident angles and polarizations, especially near the resonance frequency. This means that the materials can exhibit different internal properties (e.g. permittivity or surface impedance) when radiated with different incident angles and polarizations. A consequent problem then arises when using an effective media approach: which value for the impedance (or permittivity for bulk materials) should be applied when the materials are simultaneously radiated from different incident angles and polarizations. To incorporate the angle dependence for a particular polarization, high order boundary conditions are proposed and realized by using the method of moments for ISM [8–10]. However, the polarization dependence of approximately isotropic metasurfaces (e.g. hexagonal metasurfaces) cannot be properly handled in the high order boundary conditions in [8–10]. Actually only two-dimension, single-polarization cases are studied for strip gratings in [8]. High order boundary conditions for high impedance surfaces have also been extensively discussed for periodic structures [11, 12], but these discussions mainly focus on the analytical derivation of high order boundary conditions, without including the implementation of these conditions into full-wave numerical algorithms. Furthermore the structures discussed in [11, 12] are relatively simple because of the convenience for analytical treatments.

For microwave frequencies, miniaturized unit cells (e.g. convoluted loops [13]), are employed to significantly decrease these angle and polarization effects. But in optical frequencies, it is challenging for the nano technique to realize such fine-detail metal structures with good accuracy. Therefore, homogenization models with angle and polarization dependence are especially useful for optical metasurfaces and other possible frequency bands where techniques for miniaturized unit cells are not easily available. It is also noteworthy that in some cases the wave-vector dispersion is an advantage for metamaterial designs [14].

In this paper, we construct an ISM with both frequency and wave-vector dispersion for a more general structure (compared with the simpler structures in [11,12]), a loop-type hexagonal PMS, to simultaneously incorporate the angle and polarization dependence of the surface impedance. The most important term in our new model is the kk dyadic term that can simultaneously account for several important response features of the original hexagonal structure. The constitutive relation of this model is transformed into auxiliary differential equations (ADEs) in real space and time domain. These ADEs are then integrated into the finite-difference time-domain (FDTD) method. Finally, a finite large planar sample of the hexagonal PMS is considered and studied to test our wave-vector-dispersive model and FDTD algorithm.

2. Extracting surface impedance from S parameters of a loop-type hexagonal periodic metasurface

We consider a loop-type hexagonal periodic metasurfaces composed of infinitely thin perfect electric conductor as sketched in Fig. 1. The transmittance of this structure is computed in Fig. 2. Both s and p polarizations are considered. The incident elevation angle is swept from 0° to 75° with a 15° step, and the azimuth angle is swept (but not explicitly shown) from 0° to 360° also with 15° step. The angular frequency is normalized with respect to the resonance frequency ω₀ at normal incidence. For hexagonal periodic structure, the cross-polarization effects are zero for normal incidence [15,16] and negligible for oblique incidence. Therefore, the cross-polarization effects are neglected in this paper.

Fig. 1 The structure of the hexagonal periodic metasurface. It may also be referred as frequency selective surfaces in microwave engineering. The lattice constant is 0.25λ₀ for this structure, where λ₀ is the resonant wavelength under normal incidence.

Download Full Size | PDF

Fig. 2 The s,p-polarization transmittance (in db) of the original structure in Fig. 1 for different elevation angles of incidence. Azimuth angle ϕ = 0. The results are obtained by the full-wave commercial softwave, Computer-Simulation-Technology microwave studio.

Download Full Size | PDF

The surface impedance (1/admittance) of this infinitely thin structure is defined as a boundary condition in [6],

E_{tan}^{+} = E_{tan}^{-} \equiv E_{tan},

H_{tan}^{+} - H_{tan}^{-} = J_{s} = Y E_{tan},

where E and H are the electric and magnetic fields respectively; the subscript “tan” means the tangential components; the superscripts “+” and “−” correspond to positive and negative sides of the surface; J_s is the surface current density; Y is the admittance (1/impedance). These boundary conditions are actually the zero-thickness limit of the generalized sheet transition conditions as described in [17]. According to these boundary conditions, the reflectance and transmittance can be calculated,

t_{s} = \frac{2 cos θ}{Y / Y_{0} + 2 cos θ}, r_{s} = \frac{- Y / Y_{0}}{Y / Y_{0} + 2 cos θ},

t_{p} = \frac{2}{2 + cos θ Y / Y_{0}}, r_{p} = \frac{- cos θ Y / Y_{0}}{2 + cos θ Y / Y_{0}},

where t is transmittance; r is reflectance; subscripts “s” and “p” represent the two independent polarizations;

Y_{0} = \sqrt{μ_{0} / ∊_{0}} \approx \frac{1}{376.7 Ω}

; θ is the elevation angle of incidence. We will elaborate below that for different incident angle and polarization, the admittance Y has different values. This effect will be our main topic throughout this paper.

By use of Eq. (3) and (4), the admittance can be retrieved from the S parameters in Fig. 2. The results are shown for s and p polarizations, respectively in Fig. 3. The frequency response of the admittance of a lossless system can usually be fitted to rational functions [18] with real poles according to the Foster theorem. We employ the vector fitting algorithm [18] to extract the resonance frequency (the pole) and resonance strength (the residue) for each pair (elevation and azimuth) of incident angles. Near the resonance frequency it is enough to consider only one pair of opposite poles. Then the frequency response of Y is fitted to the following form,

Y (ω, θ, ϕ) = \frac{j F (θ, ϕ) ω}{G (θ, ϕ) - ω^{2}},

where θ, ϕ are elevation, azimuth angles of incidence respectively; F and G are positive-definite functions of θ and ϕ. F is the reciprocity of self-inductance and G is the square of resonance frequency respectively, according to effective circuit model. The angle dependence of F and G can easily be transformed into wave-vector dependence. From here on we suppose that the metasurface is in y–z plane, and then we have

k_{y} = ω_{r} sin θ cos ϕ / c, k_{z} = ω_{r} sin θ sin ϕ / c,

where

ω_{r} = \sqrt{G (θ, ϕ)}

is the resonance frequency; c is speed of light in vacuum. Now we can change the variables from angles θ and ϕ to tangential wave vectors k_y and k_z. The quantities F and G as functions of k_y and k_z are plotted in Fig. 4 for s and p polarization respectively.

Fig. 3 The frequency response of admittance Y (in Y₀) for s and p polarization and for different elevation angles of incidence. The real part of Y is not significant and not shown here, because the structure is lossless. Azimuth angle ϕ = 0.

Download Full Size | PDF

Fig. 4 The quantities G (top) and F (bottom) as functions of tangential components of wave vector k_y and k_z. F and G are in unit of Y₀ω₀ and $ω_{0}^{2}$ respectively and are related to resonance strength and frequency respectively. k_y and k_z are both in unit of ω₀/c

Download Full Size | PDF

3. Impedance sheet model with both frequency and wave-vector dispersion

Figure 4 shows that this loop-type hexagonal metasurface has almost isotropic resonance properties, except for some unimportant hexagonal features at large |k|. Therefore we can construct an isotropic model with frequency and wave-vector dispersion. The frequency dispersion is already included in Eq. (5) required by Foster theorem. The wave-vector dispersion can be extracted from Fig. 4. In order to simultaneously and isotropically incorporate the polarization and incident angle effects, we extend the quantities F and G into dyadic forms,

F = F_{1} (k^{2}) I + F_{2} (k^{2}) kk,

G = G_{1} (k^{2}) I + G_{2} (k^{2}) kk,

where I is the unity dyad; F_1,2 and G_1,2 are scalar functions in the k_y − k_z space; k is short for k_tan; k = |k|; the dyads are 2-by-2 matrixes. We can observe several important points:

(1) S polarization (perpendicular to k) and p polarization (parallel to k) are eigen vectors of the dyads F and G. It means that this model conserves s and p polarization and therefore does not have cross-polarization effects.
(2) The eigen values of the dyads F and G are F₁ and G₁ for s polarization and F₁ + F₂k² and G₁ + G₂k² for p polarization. They only depend on the k², which means that our model is truly isotropic although it has a dyadic form.
(3) The dyads F and G commute with each other and therefore the fraction in the definition of admittance Y in Eq. (5) does not have ambiguity.

Notably the kk dyadic terms properly account for the discrepancies between s and p polarizations. It is this term that distinguishes our model from the one in [8–10] which cannot properly treat the polarization properties.

Considering Fig. 4, we can expand F_1,2 and G_1,2 as Taylor series of k² and only preserve leading terms,

F_{1} = f (1 + a k^{2} + a^{'} k^{4}),

F_{2} = f (b + b^{'} k^{2}),

G_{1} = g (1 + c k^{2}),

G_{2} = g d .

where f, g, a, b, c, d are scalar constants to be determined from fitting procedure. The data in Fig. 4 should then be fitted to the model defined in Eq. (9)–(12). Before that we add some explanations. The value of f and g can directly be determined from normal incidence and can be used as the inputs of the fitting procedure. According to Fig. 4, the leading corrections for F and G are negative and positive, respectively. However, we remember that F and G must be positive-definite from previous sections. A negative k² terms in F will make the quantity negative for large enough k². Then we include a higher order corrections k⁴ for F to regularize its large-k behavior. According to Fig. 4, the k⁴ correction is not numerically significant at all, so we have much freedom in determining the magnitude of the second order corrections. One convenient choice is

a^{'} = a^{2} / 4,

b^{'} = {(a + b)}^{2} / 4 - a^{'} .

This condition ensures the positiveness of the resulted model and preserves satisfactory fitting accuracy. If F or G is not positive-definite for large k², the FDTD simulations in the following sections are quite unstable. In words, the k⁴ regularization term for F is not significant for the fitting accuracy but important for the simulation stability.

Table 1 below gives the fitting values of F and G.

Table 1. Fitting results of the parameters in F and G dyadic functions. k₀ ≡ ω₀/c.

View Table

4. Validation of the fitted model

Table 1 shows that the fitting is satisfactory because the averaged adjusted R-square is above 90%. Using the fitting table, we can interpolate or extrapolate the S parameters for any given incident angles and polarizations. The steps are as follows: firstly calculate the tangential wave vector under the given incident angles by Eq. (6); secondly determine the resonant frequency and strength (characterized by the scalar functions F and G) by Eq. (9–12); thirdly determine the corresponding admittance by Eq. (5); fourthly calculate the required S parameters for the given incident angle and polarization by Eq. (3,4).

To be convinced that the new model has better accuracy than the model without wave-vector dispersion, we test the fitted model under 85° incidence and compare the results with the full-wave simulation and the model without wave-vector dispersion. The model without wave-vector dispersion is obtained by turning off the wave-vector dependent terms in our new model. We can see from Fig. 5 that our model can give better results than the model without wave-vector dispersion. This is largely because new model can properly characterize the resonant frequency shift with incidence angles. We choose a 85° incidence angle to test our model in this paper because the resonant frequency shifts by a large amount for a large angle and the wave-vector dispersion effects are pronounced.

Fig. 5 The transmittance (in db) of the structure in Fig. 1 under 85° s-polarization (upper) and p-polarization (lower) incidence.

Download Full Size | PDF

5. Auxiliary differential equation of the new impedance sheet model

FDTD algorithm including both frequency and wave-vector dispersion is proposed for metal wire media in three dimensions in [19]. In this paper, the wave-vector-dispersive model defined by Eq. (5, 7, 8) should be transformed into time-space domain in order to be implemented in FDTD algorithm. By use of the following relation

j ω \to \frac{\partial}{\partial t}, - j k \to \nabla,

we arrive at the following ADEs,

\frac{\partial}{\partial t} P = J_{s},

\frac{\partial}{\partial t} J_{s} = FE - GP .

where the vector P is an auxiliary 2-component vector field defined on metasurface y–z plane which can be regarded as the surface polarization density; the 2 × 2 dyads F and G are defined as follows,

\begin{array}{l} F_{y y} & = & f [1 - (a + b) \partial_{y}^{2} - a \partial_{z}^{2} + (a^{'} + b^{'}) \partial_{y}^{4} \\ + (2 a^{'} + b^{'}) \partial_{y}^{2} \partial_{z}^{2} + a^{'} \partial_{z}^{4}], \end{array}

F_{y z} = F_{z y} = f [- b \partial_{y} \partial_{z} + b^{'} \partial_{y}^{3} \partial_{z} + b^{'} \partial_{y} \partial_{z}^{3}],

\begin{array}{l} F_{z z} & = & f [1 - (a + b) \partial_{z}^{2} - a \partial_{y}^{2} + (a^{'} + b^{'}) \partial_{z}^{4} \\ + (2 a^{'} + b^{'}) \partial_{y}^{2} \partial_{z}^{2} + a^{'} \partial_{y}^{4}], \end{array}

and

G_{y y} = g [1 - (c + d) \partial_{y}^{2} - c \partial_{z}^{2}],

G_{y z} = G_{z y} = - g d \partial_{y} \partial_{z},

G_{z z} = g [1 - (c + d) \partial_{z}^{2} - c \partial_{y}^{2}],

6. FDTD simulation of a finite large sample

The FDTD algorithm for the impedance sheet model of frequency selective surfaces without wave-vector dispersion has been realized in [20] but has been restricted within normal incidence. We can develop beyond normal incidence by integrating the ADE of our model in Eq. (16) and (17) into three-dimension FDTD algorithm by using the boundary condition in Eq. (1) and (2). For simplicity, we assume the impedance sheet and the tangential electric fields E_y and E_z are in the same y–z plane. E_i, P_i and J_{s_i} (i = y, z) are defined on the same points in Yee cells, respectively. The finite difference form of high-order terms are defined in the Appendix. A finite large impedance sheet (9λ₀ × 9λ₀, λ₀ = 2π/k₀) is illuminated by quasi-plane wave from a phase-modulated source of finite extent (5λ₀ × 5λ₀) very close to the impedance sheet as in Fig. 6. The incident wave has 85° incident angle and p polarization. We first consider a model without wave-vector dispersion, i.e. a = b = c = d = a′ = b′ = 0, f, g ≠ 0. The frequency of the incident wave is the same with the resonance frequency ω₀. At resonant frequency, the metasurface has zero impedance (or admittance pole) and behaves like a PEC. Therefore, the incident wave is almost totally reflected for 85° incidence except that a small portion of the wave is diffracted from the edges. If we add the wave-vector dispersion as shown in Table 1, the result is quite different as shown in Fig. 6. The transmitted wave along the incident 85° direction is quite obvious because the resonance frequency is already shifted for 85° incidence. As references, the actual transmittance rate (in db) of the original structure for 85° incident angle and ω₀ frequency are about −8db for p polarization according to Fig. 5. Figure 7 shows the s polarization case under the same conditions as Fig. 6. We may find in Fig. 7 that the models with and without wave-vector dispersion have small difference. According to Fig. 5 the s-polarization transmittance for the real structure is −40db, which is consistent with the results in Fig. 7. Therefore, our new FDTD algorithm succeeds in incorporating the angle and polarization dependences of the PMS simultaneously.

Fig. 6 The upper (lower) figure is the case without (with) wave-vector dispersion. The incident wave has p polarization. Perfect matching layers are used. The space discretization step is δx = δy = δz = 0.1λ₀. Time discretization step is 0.6 times the Courant limit for free vacuum, δt = 0.0346λ₀/c. The z components of the magnetic fields are shown in the two figures. The black solid line represents the position of the source.

Download Full Size | PDF

Fig. 7 The upper (lower) figure is the case without (with) wave-vector dispersion. The incident wave has s polarization. Perfect matching layers are used. The space discretization step is δx = δy = δz = 0.1λ₀. Time discretization step is 0.6 times the Courant limit for free vacuum, δt = 0.0346λ₀/c. The z components of the electric fields are shown in the two figures. The black solid line represents the position of the source.

Download Full Size | PDF

To evaluate the accuracy of the diffracted wave from edges in our model compared with the cases in the real structure is a complicated problem, as we must perform a full-wave study of a large matrix of unit cells which is really computationally expensive. However, our model does give a better characterizing for the transmitted wave in the direction of incidence in Fig. 6. For a even larger metasurface sample which is much larger than the lateral width of the incident beam, the edge diffraction may be neglected and our model can work better.

7. Conclusions

We propose a generalized impedance sheet model with both frequency and wave-vector dispersions for a loop-type periodic hexagonal metallic metasurfaces that have approximately isotropic responses near resonance. This model can properly and accurately considers the incident angle and polarization dependence of such metasurface, especially the shift of resonance frequency for different incident angles and polarizations. The integration of this model into FDTD algorithm with the aid of auxiliary differential equations enables us to numerically analyze all the angle and polarization effects in a single run. The integration of our model into method of moments or finite elements is to be studied. Our model can be generalized to study other types of periodic metasurfaces with different symmetries, (e.g. anisotropic structures with cross-polarization effects).

8. Appendix: finite difference expressions for the high-order terms

In the following, $E_{y}^{(z z)}$ is short for $\partial_{z}^{2} E_{y}$ . For other cases, the notations are similar.

In our implementation, the field second-order derivatives $E_{y}^{(z z)}$ , $E_{y}^{(y y)}$ , $E_{y}^{(y z)}$ are defined on the same points as E_y; $E_{z}^{(z z)}$ , $E_{z}^{(y y)}$ , $E_{y}^{(y z)}$ are defined on the same points as E_z.

The typical second-order terms are,

\begin{array}{r} E_{y}^{(z z)} (n_{y} + 1 / 2, n_{z}) = [E_{y} (n_{y} + 1 / 2, n_{z} + 1) + \\ E_{y} (n_{y} + 1 / 2, n_{z} - 1) - 2 E_{y} (n_{y} + 1 / 2, n_{z})] / d z^{2} . \end{array}

\begin{array}{r} E_{y}^{(y z)} (n_{y} + 1 / 2, n_{z}) = [E_{z} (n_{y} + 1, n_{z} + 1 / 2) + \\ E_{z} (n_{y}, n_{z} - 1 / 2) - E_{z} (n_{y} + 1, n_{z} - 1 / 2) \\ - E_{z} (n_{y}, n_{z} + 1 / 2)] / d z^{2} . \end{array}

The fourth-order terms can be defined as the second-order differences of the second-order terms, e.g.

E_{z}^{(y z z z z)} (n_{y} + 1 / 2, n_{z}) = E_{z}^{(y z) (z z)} .

Similar definitions can be applied to the P field in Eq. (16) and (17).

Funding

Yantai University (No. WL17B01); National Key Research and Development Program of China (No. 2016YFB0402705); and Shenzhen Key Lab Fund (No. ZDSYS20170228105421966).

References and links

1. A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, “Planar photonics with metasurfaces,” Science 339(15), 1232009 (2013). [CrossRef] [PubMed]

2. B. A. Munk, Frequency Selective Surfaces: Theory and Design (John Wiley Sons, Inc., 2000). [CrossRef]

3. Y. Zhao, X.-X. Liu, and A. Alu, “Recent advances on optical metasurfaces,” Journal of Optics 16, 123001 (2014). [CrossRef]

4. R. Mittra, “A look at some challenging problems in computational electromagnetics,” IEEE Antennas and Propagation Magazine 46(5), 18–32 (2004). [CrossRef]

5. Y. Rahmat-Samii and A. N. Tulintseff, “Diffraction analysis of frequency selective reflector antennas,” IEEE Trans. Antennas Propag. 41(4), 476–487 (1993). [CrossRef]

6. K. W. Whites and R. Mittra, “An equivalent boundary-condition model for lossy planar periodic structures at low frequency,” IEEE Trans. Antennas Propag. 44(12), 1617–1629 (1996). [CrossRef]

7. Y. Zhao, N. Engheta, and A. Alu, “Homogenization of plasmonic metasurfaces modeled as transmission-line loads,” Metamaterials , 5(2), 90–96 (2011). [CrossRef]

8. B. Stupfel and Y. Pion, “Impedance boundary conditions for finite planar or curved frequency selective surfaces,” IEEE Trans. Antennas Propag. 53(4), 1415–1425 (2005). [CrossRef]

9. B. Stupfel, “Impedance boundary conditions for finite planar or curved frequency selective surfaces embedded in dielectric layers,” IEEE Trans. Antennas Propag. 53(11), 3654–3663 (2005). [CrossRef]

10. B. Stupfel, “Implementation of high-order impedance boundary conditions in some integral equation formulations,” IEEE Trans. Antennas Propag. 63(4), 1658–1668 (2015). [CrossRef]

11. O. Luukkonen, G. Granet, G. Goussetis, D. Lioubtchenko, A. V. Raisanen, and S. A. Tretyakov, “Simple and accurate analytical model of planar grids and high-impedance surfaces comprising metal strips or patches,” IEEE Trans. Antennas Propag. 56(6), 1624–1632 (2008). [CrossRef]

12. S. A. Tretyakov, Analytical Modeling in Applied Electromagnetics (Artech House, 2003).

13. F. Huang, J. Batchelor, and E. Parker, “Interwoven convoluted element frequency selective services with wide bandwidths,” Electronics Letters 42(14), 788–790 (2006). [CrossRef]

14. J. Luo, Y. Yang, Z. Yao, W. Lu, B. Hou, Z. H. Hang, C. Chan, and Y. Lai, “Ultratransparent media and transformation optics with shifted spatial dispersions,” Phys. Rev. Lett. 117, 223901 (2016). [CrossRef] [PubMed]

15. A. Mackay, “Proof of polarisation independence and nonexistence of crosspolar terms for targets presenting n-fold (n>2) rotational symmetry with special reference to frequency-selective surfaces,” Electron. Lett. 25(24), 1624–1625, (1989). [CrossRef]

16. K. W. Whites and R. Mittra, “Scattered field properties of symmetric periodic structures,” IEEE Trans. Antennas Propag. 42(5), 722–727 (1994). [CrossRef]

17. C. L. Holloway, E. F. Kuester, J. A. Gordon, J. O’Hara, J. Booth, and D. R. Smith, “An overview of the theory and applications of metasurfaces: The two-dimensional equivalents of metamaterials,” IEEE Antennas and Propagation Magazine 54(2), 10–35 (2012). [CrossRef]

18. B. Gustavsen and A. Semlyen, “Rational approximation of frequency domain responses by vector fitting,” IEEE Transactions on Power Delivery 14(3), 1052–1061 (1999). [CrossRef]

19. Y. Zhao, P. A. Belov, and Y. Hao, “Modelling of wave propagation in wire media using spatially dispersive finite-difference time-domain method: numerical aspects,” IEEE Trans. Antennas Propag. 55(6), 1506–1513 (2006). [CrossRef]

20. M. K. Karkkainen and P. M. T. Ikonen, “Finite-difference time-domain modeling of frequency selective surfaces using impedance sheet conditions,” IEEE Trans. Antennas Propag. 53(9), 2928–2937 (2005). [CrossRef]

Item	Fitting value	Adjusted R-square
f/(Y₀ω₀)	1.0	—
$a k_{0}^{2}$	−0.10	0.99
$b k_{0}^{2}$	0.066	0.92
$g / ω_{0}^{2}$	1	—
$c k_{0}^{2}$	0.043	0.99
$d k_{0}^{2}$	−0.022	0.90

Wave-vector and polarization dependent impedance model for a hexagonal periodic metasurface exemplified through finite-difference time-domain simulations

Abstract

1. Introduction

2. Extracting surface impedance from S parameters of a loop-type hexagonal periodic metasurface

3. Impedance sheet model with both frequency and wave-vector dispersion

4. Validation of the fitted model

5. Auxiliary differential equation of the new impedance sheet model

6. FDTD simulation of a finite large sample

7. Conclusions

8. Appendix: finite difference expressions for the high-order terms

Funding

References and links

Cited By

Figures (7)

Tables (1)

Equations (26)

Optics Express