Ultrashort elliptically polarized laser pulse interaction with helical photonic metamaterial

N. N. Potravkin; E. B. Cherepetskaya; I.A. Perezhogin; V.A. Makarov

doi:10.1364/OME.4.002090

1. Introduction

During last years metamaterials representing themselves a 2D-lattice composed of 3D-helices attract a great interest owing to their application in compact devices for polarization control of ultrashort laser pulses (see [1] and references therein). Such structures have a big difference between the transmission coefficients for circularly polarized light pulses with opposite handedness of polarization rotation. They are very prospective [1] for the generation of circularly polarized light, which required in various [2–4] applications. The traditional methods of generation of monochromatic circularly polarized light, based on the utilization of a $λ / 4 -$ plate (birefringent plate with a thickness equal to one fourth of a wavelength) or the cholesteric liquid crystal with pitch of its helical structure close to the wavelength of the propagating light cannot be used for the ultrashort pulses with broad spectrum.

In [1] it was experimentally demonstrated that a metamaterial consisting of metallic 3D helices can be used as a broad-band thin-film polarizer for electromagnetic radiation, transforming the polarization of the transmitted light into circular. Spectral properties of this object were analyzed using the commercial calculating software [5, 6]. The unit cell of such a material can be composed of two coaxial helices oriented in a peculiar way one relatively to the other. This allows one to make the operating frequency range of such a polarizing device broader at almost 50% [6] comparing to conventional metamaterial polarizer [1]. If the unit cell contains bigger amount of the coaxial helices, then the signal-to-noise ratio significantly improves in such a device [5]. In [7] the polymeric material with 8-periods helices was used. It was shown that this sample provides twenty times difference between the mean transmission coefficients for the circularly polarized light with the opposite handedness of the polarization rotation, if the angle of incidence (relatively to the axes of the helices) was less than 7 degrees.

However, in the abovementioned works the authors consider mainly the spectral features of the transmitted and reflected radiation, while the polarization evolution of the propagating pulses (in case of relatively long ones) and the temporal behavior of the electric field vector in them (in ultrashort pulses) remained outside their scope. In addition, for the polymeric metamaterials (unlike for the metallic nano- and microstructures), the light-matter interaction and the role of the geometry of the metamaterial unit cell was not studied so thoroughly.

In the present work the finite-difference time-domain (FDTD) method is applied to the study of polymeric metamaterial unit cell analogous to that in [7], where the effect of polarization-selective transmission was observed experimentally. We model the problem of light propagation in such a metamaterial numerically, and particularly we study the influence of metamaterial unit cell parameters (first of all, of a number of helix coils) on the transmission and reflection of the elliptically polarized light. In case of long pulses with duration of about 20 light-field oscillations or more the temporal dynamics of polarization of transmitted and reflected light can be described analogously to [8, 9] by the hodograph of the electric field vector in the electromagnetic wave. In case of ultrashort pulses with much broader spectra we define the parameters characterizing the polarization in spectral domain, because the interpretation of the hodograph of the electric field becomes too complex.

We discuss the origins of the optical resonances causing the strong difference in optical properties of the metamaterial under consideration for the transmission of the RHCP and LHCP light. Unlike in metallic 3D helix metamaterials, where the optical resonance is caused by the oscillations of the conduction electrons in metal [5], the polymeric dielectric material causes the specific oscillations of the electric and magnetic parts of the electromagnetic field energy.

2. Formulation of the problem

Let us consider the plane electromagnetic wave propagating along the $z -$ axis in a thin layer of a periodic medium. This medium periodic in x- and y-directions can be represented by a unit cell. Each cell contains a right-handed (or left-handed) helix with $n$ coils having pitch $h$ (Fig. 1).Such a helix can be manufactured from the isotropic dielectric material with dielectric permittivity $ε_{1}$ using a technique described in [5, 6]. The $z -$ axis is parallel to the axis of the helix.

Fig. 1 Schematic image of a helix, representing the metamaterial unit cell. The definition of characteristic sizes is shown in the scheme.

Download Full Size | PDF

Maxwell equations and material equations connecting strength and induction of electric ( $E (x, y, z, t)$ and $D (x, y, z, t)$ ) and magnetic ( $H (x, y, z, t)$ and $B (x, y, z, t)$ ) fields can be written as following:

\frac{1}{c} \frac{\partial B_{x}}{\partial t} = \frac{\partial E_{y}}{\partial z} - \frac{\partial E_{z}}{\partial y} \frac{1}{c} \frac{\partial D_{x}}{\partial t} = \frac{\partial H_{z}}{\partial y} - \frac{\partial H_{y}}{\partial z} \frac{1}{c} \frac{\partial B_{y}}{\partial t} = \frac{\partial E_{z}}{\partial x} - \frac{\partial E_{x}}{\partial z},

\frac{1}{c} \frac{\partial D_{y}}{\partial t} = \frac{\partial H_{x}}{\partial z} - \frac{\partial H_{z}}{\partial x} \frac{1}{c} \frac{\partial B_{z}}{\partial t} = \frac{\partial E_{x}}{\partial y} - \frac{\partial E_{y}}{\partial x} \frac{1}{c} \frac{\partial D_{z}}{\partial t} = \frac{\partial H_{y}}{\partial x} - \frac{\partial H_{x}}{\partial y},

D_{i} = ε (x, y, z) E_{i}, B_{i} = H_{i}

In Eq. (3) $ε (x, y, z) = ε_{1}$ , if the point with coordinates $(x, y, z)$ lie within the helix area and $ε (x, y, z) = 1$ if $(x, y, z)$ is in empty space. We consider the elliptically polarized light pulse with central frequency of its spectrum $ω_{0}$ , half width of its profile (along z-axis) w₀ and the wavelength $λ$ . Initially ( $t = 0$ ) the light pulse is incident on the border of the medium located at $z = - n h / 2$ . The Cartesian components of the electric field in this light pulse are given by the following expressions:

\begin{array}{l} E_{x} (t = 0) = (^{I_{0} / 2) 1 / 2} [^{1 - (^{1 - M_{0}^{2}) 1 / 2}] 1 / 2} \exp [- (^{z - z_{0}) 2} / w_{0}^{2}] \times \\ \times s i g n {M_{0}} \sin (2 π (z - z_{0}) / λ), \end{array}

\begin{array}{l} E_{y} (t = 0) = (^{I_{0} / 2) 1 / 2} [^{1 + (^{1 - M_{0}^{2}) 1 / 2}] 1 / 2} \exp [- (^{z - z_{0}) 2} / w_{0}^{2}] \times \\ \times \cos (2 π (z - z_{0}) / λ) . \end{array}

At $z = z_{0}$ the maximum value of their dimensionless intensity $I (z) = [E_{x}^{2} (z) + E_{y}^{2} (z)] / I_{0}$ is achieved. The value of the ellipticity degree of the polarization ellipse $M_{0}$ of the incident pulse can be chosen from $1$ to $- 1$ (see also [10]). $M_{0} = 1$ (or $M_{0} = - 1$ ) corresponds to the right- handed (or left-handed) circular polarization, and $M_{0} = 0$ corresponds to the linear polarization.

In our calculations we used the value of time step as $0.043$ fs and the value of space step $26$ nm. As in [7] we considered the translation period (in $x -$ and $y -$ directions) of the metamaterial unit cell (the lattice constant of the metamaterial) as $1.3$ μm, the helix pitch as $1.3$ μm, the spiral diameter as $0.79$ μm, the lateral diameter of the spiral arms as $0.38$ μm, the axis diameter of the spiral arms as $0.83$ μm. All the above mentioned parameters are shown in the scheme in Fig. 1. The dielectric permittivity of the helix was taken as $ε = 2.47$ . The numerical methods, their substantiation and the algorithm of the numerical scheme are described in the Appendix.

3. The discussion of results

In case of a long pulse ( $w = 20 λ$ , $λ = 1.67$ μm, $z_{0} = - 50 λ$ ), and if the helices have big number of coils n there takes place selective reflection of the circularly polarized light. Namely, in case of right-handed helices, the LHCP light passes through the metamaterial, while the RHCP light appears to be almost wholly reflected. And we observe the opposite behavior for the left-handed helices. In both cases the z-components of the electric field vector in the reflected and transmitted waves turn to zero for $z < < - n h / 2$ and for $z > > n h / 2$ $z > > n h / 2$ (at long distance from the sample), while $E_{x, y} = E_{x, y} (x, y, z)$ components of these waves practically do not change in the $x y -$ plane. In Fig. 2 there are shown typical hodographs of the electric field vector (trajectory traced by the end of vector $E (z)$ in $x y z$ -space) in light pulses reflected from (a, c and e) and transmitted through (b, d and f) the metamaterial for different values of $M_{0}$ in case of right-handed helices. Designations ${E^{'}}_{x, y} = E_{x, y} / I_{0}^{1 / 2}$ are used in Fig. 2.

Fig. 2 Hodographs of the electric field vector in light pulses reflected from (a, c and e) and transmitted through (b, d, and f) the metamaterial, $t = 1075$ fs, $n = 8$ , $M_{0} = 0$ (а, b), $M_{0} = 1$ (c, d), $M_{0} = - 1$ (e, f).

Download Full Size | PDF

In case of the incidence of the linearly polarized pulse the peak intensities of the reflected and transmitted pulses are approximately the same (Figs. 2(a) and 2(b)). In case of the incidence of the RHCP light pulse ( $M_{0} = 1$ ) onto a medium composed of right-handed helices (Figs. 2(c) and 2(d)), the peak value of the electric field strength in the transmitted pulse is order of magnitude smaller than that in case of the incidence of LHCP pulse (Figs. 2(e) and 2(f)). In the latter case the transmitted radiation is almost circularly polarized (lef-hand), and the elliptically polarized reflected pulse has complicated profile (Fig. 2(e)).

In order to describe the evolution of polarization of long pulse we use an approach formulated in [8, 9]. Within this approach instead of continuous dependencies $E_{x, y} (z)$ we use a sequences of polarization ellipses determined by the values $E_{x, y}$ in points $z = {\tilde{z}}_{m}$ , where the intensity I achieves local maxima. The ellipticity degree and the angle of orientation of the polarization ellipses in these points can be found using the following formulae:

| M ({\tilde{z}}_{m}) | = \frac{2^{1 / 2} I^{1 / 2} ({\tilde{z}}_{m}) {[I ({\bar{z}}_{m}) + I ({\bar{z}}_{m + 1})]}^{1 / 2}}{I ({\tilde{z}}_{m}) + [I ({\tilde{z}}_{m}) + I ({\tilde{z}}_{m + 1})] / 2},

Ψ ({\tilde{z}}_{m}) = - a r c t g [E_{x} ({\tilde{z}}_{m}) / E_{y} ({\tilde{z}}_{m})],

Here

{\bar{z}}_{m}

are the points of local minima of

I (z)

, their numbering m is given in such a way that

{\bar{z}}_{m} \leq {\tilde{z}}_{m} \leq {\bar{z}}_{m + 1}

. The sign of

M ({\tilde{z}}_{m})

is determined by the handedness of rotation of the electric field vector. Defined above discrete dependencies

M ({\tilde{z}}_{m})

and

Ψ ({\tilde{z}}_{m})

show the evolution of the polarization ellipse ellipticity degree and orientation angle along the pulse envelope. In the extreme case, when the monochromatic radiation propagates, the defined above quantities coincide with the ellipticity degree and the angle of orientation of the polarization ellipse of monochromatic light.

In case of incidence of linearly polarized pulse, the transmitted pulse has bell-like shape. If we increase the number $n$ of helix coils (Fig. 3), the intensity peak in the pulse shifts and the maximum intensity value becomes almost 30% less. If $n = 2$ the transmitted pulse is elliptically polarized, and if $n = 8$ the transmitted pulse is almost LHCP (see the black and the blue curves in Fig. 3(a)). In latter case the z-dependence of the angle of orientation of the polarization ellipse (along the pulse envelope) becomes non-monotone (Fig. 3(b)). When the RHCP pulse falls onto a medium consisting of right-handed helices, the intensity of the transmitted radiation decreases exponentially with the increase of n, and its polarization becomes closer to LHCP (Fig. 3(c)). In addition, the monotone change of orientation angle along the pulse becomes steeper with the increase of $n$ (Fig. 3(d)). The $π -$ jump of $Ψ -$ dependence is due to the definition of this quantity (see Eq. (5)). In case of incidence of LHCP pulse, the peak intensity is the biggest and the ellipticity degree of the polarization ellipse is most close to the initial one in case of $n = 2$ (Fig. 3(e)). In this case the increase of the number of helical coils results in non-monotone character of the change of polarization ellipse orientation angle along the pulse (Fig. 3(f)).

Fig. 3 The dependencies of the ellipticity degree $M$ (a, c, and e) and the angle of orientation of the polarization ellipse $Ψ$ (b, d and f) on the longitudinal coordinate z along the temporal envelope of a pulse traversed through the metamaterial consisting of right-handed helices. The dependencies are shown at time moment $t = 1075$ fs (in a far-field region after the metamaterial) in case of the incidence of (a and b) – linearly polarized pulse; (c and d) – RHCP pulse; (e and f) – LHCP pulse. Black curves correspond to $n = 2$ ; red curves: $n = 4$ ; blue curves: $n = 8$ .

Download Full Size | PDF

Our studies have shown that in case of incidence of long circularly polarized pulse on a border of the metamaterial, there can take place various modes of oscillation of electric and magnetic constituents of the energy of the electromagnetic field. Let us consider the RHCP pulse incidence on the metamaterial consisting of right-hand helices. In this case there are such time values $t_{p}$ ( $p = 1, 2, 3, \dots$ ) separated with a period of 2π/ω, when the density of the electric field energy $w_{e} (t, x, y, z) = (D \cdot E) / 8 π$ becomes concentrated (non-zero) only in the helices (Fig. 4(a)), while the density of the magnetic field energy $w_{h} (t, x, y, z) = (B \cdot H) / 8 π$ turns to zero within the whole metamaterial (both within and outside the helices). Simultaneously $w_{e}$ exponentially decays as the field propagates further in a bulk of a metamaterial (Figs. 4(a) and 4(b)). After half of the period passes, $w_{e} = 0$ in a whole metamaterial, while the magnetic part of the energy is concentrated outside the helices between them (Fig. 4(b)), and also decays exponentially in the course of propagation into the metamaterial.

Fig. 4 Spatial distributions of $w_{e} (x, y = 0, z)$ and $w_{h} (x, y = 0, z)$ in case of propagation of the long pulses at time values $t_{p}$ (а), (c) and after half of the oscillation period, at $t_{p} + π / ω$ (b), (d). RHCP incidence – (a) and (b); LHCP incidence – (c) and (d). The metamaterial consists of right-hand helices, $n = 8$ .

Download Full Size | PDF

When the LHCP pulse falls on the metamaterial consisting of right-handed helices the character of the oscillations of the electric and magnetic field energy is essentially different from the described above for the RHCP pulse. The values of $w_{e}$ and $w_{h}$ are always non-zero in any point of the metamaterial, and after a half of the oscillation period both electric and magnetic field energies move from one end of the helix to another (Figs. 4(c) and 4(d)). Unlike in previous case the magnetic part of the energy now is concentrated in a space between the helices.

If one decreases the duration of the incident pulse the abovementioned processes make the interpretation of hodographs of the electric field in transmitted and reflected pulses very complex and almost unsuccessful (Fig. 5).The discrete dependencies $M ({\tilde{z}}_{m})$ and $Ψ ({\tilde{z}}_{m})$ defined in (6), (7) also are changing irregularly with the quasi-period of oscillations comparable with the wavelength. The information about the polarization of the transmitted radiation can be accessed through the transmission coefficient $T (ω) = {(| S_{x}^{t} |^{2} + | S_{y}^{t} |^{2})}^{1 / 2} \cdot {(| S_{x}^{i} |^{2} + | S_{y}^{i} |^{2})}^{- 1 / 2}$ for the radiation at given frequency ω and the ellipticity degree of its polarization ellipse $M (ω) = i (S_{y}^{t} S_{x}^{t *} - S_{x}^{t} S_{y}^{t *}) \cdot (| S_{x}^{t} |^{2} + | S_{y}^{t} |^{2})$ , which changes in a range from $- 1$ to $1$ . Here $S_{x, y}^{i} (ω)$ and $S_{x, y}^{t} (ω)$ are the Fourier images of the temporal dependencies of Cartesian components of the electric field vector in transmitted and reflected pulses (the asterisk in the superscript stands for the complex conjugated quantity).

Fig. 5 The hodographs of the electric field vectors in (a) reflected and (b) transmitted pulses at $t = 1075$ fs after the incidence of ultrashort ( $w_{0} = 2 λ$ ) linearly polarized pulse on the metamaterial; $z_{0} = - 12 λ$ .

Download Full Size | PDF

Figure 6 shows the dependencies $T (ω / ω_{0})$ and $M (ω / ω_{0})$ , where $ω_{0} = 1.16 \cdot 10^{15}$ rad/s, in case of the incidence of RHCP (red curves) or LHCP (blue curves) pulse, when a metamaterial consists of right-handed helices with $n = 8$ (solid lines) or $n = 4$ (dashed lines). If the incident radiation is RHCP, then there is a relatively broad frequency range around $ω \approx ω_{0}$ where the intensity of the transmitted RHCP radiation is only a few percents of the intensity of the incident RHCP radiation (Fig. 6(a)) and $M (ω / ω_{0}) \approx - 1$ (Fig. 6(b)). With the increase of the number of coils in the helices this frequency range becomes broader. The LHCP radiation within this frequency range almost does not changes its polarization after propagation in a metamaterial.

Fig. 6 (a) The frequency spectrum of the transmission coefficient and (b) the frequency spectrum of the ellipticity degree of the polarization ellipse in case of the incidence of RHCP (red curves) or LHCP (blue curves) pulse, when a metamaterial consists of right-handed helices with $n = 8$ (solid lines) or $n = 4$ (dashed lines). $ω_{0} = 1.16 \cdot 10^{15}$ rad/s.

Download Full Size | PDF

However, for higher frequencies it becomes linearly polarized and then with the further increase of the frequency the ellipticity degree of the polarization ellipse in spectrum domain non-monotonically decreases to $- 1$ (Fig. 6(b)). Slight and moderate variations of the following parameters: the dielectric permittivity of the helix material, the translation period of the unit cell of a metamaterial, pitch and diameter of a helix, lateral and axial diameters of the helix arms, do not change qualitatively the dependencies shown in Fig. 6, and only shift the central frequency of the “anomalous transmission” range and change the size of this range. For example, with the increase of $ε$ (when the other parameters are fixed) this range broadens and its center shifts to lower frequencies.

4. Conclusion

Using the FDTD method we have studied the influence of the parameters of polymeric metamaterial unit cell on the transmission and reflection of normally incident elliptically polarized light. In case of long incident pulses (with duration of about $20$ oscillation periods or more) we analyzed hodographs of the electric field vectors in order to describe the temporal evolution of the transmitted and reflected light pulses. If the polarization of the incident pulse is close to linear, then the intensities of the transmitted and reflected elliptically polarized pulses are approximately equal. In case of the incidence of the RHCP pulse onto a metamaterial consisting of right-handed helices, the peak intensity in the transmitted pulse is order of magnitude lower than in case of the incidence of the LHCP pulse. If the polarization of the incident radiation is close to left-hand circular polarization, than after propagation through the metamaterial the polarization of this pulse remains practically unchanged. At the same time, the reflected pulse which is elliptically polarized has rather complicated shape.

It is shown that in case of incidence of long circularly polarized pulse on a border of the metamaterial, there can take place various modes of oscillation of electric and magnetic constituents of the energy of the electromagnetic field. In case of incidence of the RHCP pulse there are such time frames, when the density of the electric field energy has non-zero value only within the helices, while the density of the magnetic field energy turns to zero within whole metamaterial (both within and outside the helices). After half of the period passes, the density of the electric field energy turns to zero in whole metamaterial, while the magnetic part of the energy becomes concentrated outside the helices between the coils. In case of incidence of the LHCP radiation on the metamaterial consisting of right-handed helices the electric and magnetic parts of the energy are always non-zero in any point of the metamaterial, and after a half of the oscillation period both electric and magnetic field energies move from one end of the helix to another.

If the incident ultrashort pulse is RHCP and the metamaterial contains of right-handed helices, then there is a reasonably broad frequency range, where the intensity of the transmitted radiation is only a few percent of that of the incident, and its polarization is close to LHCP. With the increase of the number of helix coils the size of this frequency range grows. If the LHCP radiation is incident on the metamaterial it practically does not changes its polarization in the course of propagation, however, at higher frequencies in spectrum domain its polarization becomes linear, and at even higher frequencies it becomes circular. Slight and moderate variations of the dielectric permittivity of the helix material, the translation period of the unit cell of a metamaterial, pitch and diameter of a helix, lateral and axial diameters of the helix arms, shift the central frequency of the “anomalous transmission” range and change the size of this range. However, they do not change qualitatively the overall character of light polarization evolution.

5. Appendix

In order to perform discrete representations of the spatial derivative operators of the Maxwell equations we used FDTD method, which implies that the Cartesian components of the strength of the electric and magnetic fields are calculated in the centers of the cells of the grid (Liu scheme [11]). Unlike in more widespread approach, when the Cartesian components of these vectors are calculated at different grids shifted one relatively to another, our approach is more prospective for its further implementation in nonlinear optics problems, for example, in propagation of the ultrashort laser pulse through the nonlinear metamaterial. Unlike in a linear problem, in this case the material equations do not split to independent equations for each Cartesian component of the field. Moreover, if the Cartesian components of vectors of the electric and magnetic fields are calculated not in the same point of the cell, then additional interpolation procedures are needed in order to resolve the material equations, which inevitably result in losses in precision and efficiency of the numerical scheme. As in scheme in [11] we used non-symmetric expressions for the spatial derivatives at unstaggered grid. It allows one to get rid of numerical oscillations due to odd–even decoupling [12], appearing when using symmetric (central) difference approximations in such grids. The increase of the order of the approximation allow us to overcome typical disadvantages of second-order approximation schemes, such as large numerical dispersion and scheme anisotropy on the unstaggered grid.

Let us consider homogeneous Cartesian grid and a template consisting of $r + l + 1$ points on it ( $l \geq r$ ). In order to calculate the spatial derivatives $H_{α}$ we used upwind-biased differencing operator [11]:

\nabla_{β} H_{α} (μ Δ x, ν Δ y, ρ Δ z, (σ + 1 / 2) Δ t) = \nabla_{_{β}} H_{α, μ, ν, ρ}^{(σ + 1 / 2)} = (1 / Δ β) \sum_{m = - l}^{r} a_{m} H_{α, μ + m δ_{μ β}, ν + m δ_{ν, ρ} + m δ_{ρ β}}^{(σ + 1 / 2)} .

Here α and β take values

x, y, z

. Indices

μ

,

ν

,

ρ

and

σ

determine the number of the calculation cell in

x y z t -

space,

Δ β

is a step at

β -

axis,

Δ t

is a step at time axis. In order to calculate spatial derivatives

E_{α}

we used downwind-biased differencing [11]:

Δ_{β} E_{α} (μ Δ x, ν Δ y, ρ Δ z, σ Δ t) = Δ_{β} E_{α, μ, ν, ρ}^{(σ)} = - (1 / Δ β) \sum_{m = - r}^{l} a_{- m} E_{α, μ + m δ_{μ β}, ν + m δ_{ν β}, ρ + m δ_{ρ β}}^{(σ)} .

Here coefficients

a_{m}

in (8) and (9) correspondingly must satisfy the following systems of equations:

\sum_{m = - l}^{r} a_{m} = 0

,

\sum_{m = - l}^{r} m a_{m} = 1

, …,

\sum_{m = - l}^{r} m^{n} a_{m} = 0

, where

n = 2, 3, \dots, r + l

. For the scheme with third-order approximation on the spatial coordinates we will obtain the following difference analogues of the Maxwell Eqs. (1), (2):

\begin{array}{l} \frac{H_{x, μ, ν, ρ}^{(σ + 1 / 2)} - H_{x, μ, ν, ρ}^{(σ - 1 / 2)}}{c Δ t} = - \frac{a_{- 2} E_{y, μ, ν, ρ + 2}^{(σ)} + a_{- 1} E_{y, μ, ν, ρ + 1}^{(σ)} + a_{0} E_{y, μ, ν, ρ}^{(σ)} + a_{1} E_{y, μ, ν, ρ - 1}^{(σ)}}{Δ z} + \\ + \frac{a_{- 2} E_{z, μ, ν + 2, ρ}^{(σ)} + a_{- 1} E_{z, μ, ν + 1, ρ}^{(σ)} + a_{0} E_{z, μ, ν, ρ}^{(σ)} + a_{1} E_{y, μ, ν - 1, ρ}^{(σ)}}{Δ y}, \end{array}

\begin{array}{l} \frac{H_{y, μ, ν, ρ}^{(σ + 1 / 2)} - H_{y, μ, ν, ρ}^{(σ - 1 / 2)}}{c Δ t} = - \frac{a_{- 2} E_{z, μ + 2, ν, ρ}^{(σ)} + a_{- 1} E_{z, μ + 1, ν, ρ}^{(σ)} + a_{0} E_{z, μ, ν, ρ}^{(σ)} + a_{1} E_{z, μ - 1, ν, ρ}^{(σ)}}{Δ x} + \\ + \frac{a_{- 2} E_{x, μ, ν, ρ + 2}^{(σ)} + a_{- 1} E_{x, μ, ν, ρ + 1}^{(σ)} + a_{0} E_{x, μ, ν, ρ}^{(σ)} + a_{1} E_{x, μ, ν, ρ - 1}^{(σ)}}{Δ z}, \end{array}

\begin{array}{l} \frac{H_{z, μ, ν, ρ}^{(σ + 1 / 2)} - H_{z, μ, ν, ρ}^{(σ - 1 / 2)}}{c Δ t} = - \frac{a_{- 2} E_{x, μ, ν + 2, ρ}^{(σ)} + a_{- 1} E_{x, μ, ν + 1, ρ}^{(σ)} + a_{0} E_{x, μ, ν, ρ}^{(σ)} + a_{1} E_{x, μ, ν - 1, ρ}^{(σ)}}{Δ y} + \\ + \frac{a_{- 2} E_{y, μ + 2, ν, ρ}^{(σ)} + a_{- 1} E_{y, μ + 1, ν, ρ}^{(σ)} + a_{0} E_{y, μ, ν, ρ}^{(σ)} + a_{1} E_{y, μ - 1, ν, ρ}^{(σ)}}{Δ x}, \end{array}

\begin{array}{l} \frac{D_{x, μ, ν, ρ}^{(σ + 1)} - D_{x, μ, ν, ρ}^{(σ)}}{c Δ t} = \frac{a_{- 2} H_{z, μ, ν - 2, ρ}^{(σ + 1 / 2)} + a_{- 1} H z_{z, μ, ν - 1, ρ}^{(σ + 1 / 2)} + a_{0} H_{z, μ, ν, ρ}^{(σ + 1 / 2)} + a_{1} H_{z, μ, ν + 1, ρ}^{(σ + 1 / 2)}}{Δ y} - \\ - \frac{a_{- 2} H_{y, μ, ν, ρ - 2}^{(σ + 1 / 2)} + a_{- 1} H_{y, μ, ν, ρ - 1}^{(σ + 1 / 2)} + a_{0} H_{y, μ, ν, ρ}^{(σ + 1 / 2)} + a_{1} H_{y, μ, ν, ρ + 1}^{(σ + 1 / 2)}}{Δ z}, \end{array}

\begin{array}{l} \frac{D_{y, μ, ν, ρ}^{(σ + 1)} - D_{y, μ, ν, ρ}^{(σ)}}{c Δ t} = \frac{a_{- 2} H_{x, μ, ν, ρ - 2}^{(σ + 1 / 2)} + a_{- 1} H z_{x, μ, ν, ρ - 1}^{(σ + 1 / 2)} + a_{0} H_{x, μ, ν, ρ}^{(σ + 1 / 2)} + a_{1} H_{x, μ, ν, ρ + 1}^{(σ + 1 / 2)}}{Δ z} - \\ - \frac{a_{- 2} H_{z, μ - 2, ν, ρ}^{(σ + 1 / 2)} + a_{- 1} H_{z, μ - 1, ν, ρ}^{(σ + 1 / 2)} + a_{0} H_{z, μ, ν, ρ}^{(σ + 1 / 2)} + a_{1} H_{z, μ + 1, ν, ρ}^{(σ + 1 / 2)}}{Δ x}, \end{array}

\begin{array}{l} \frac{D_{z,, μ, ν, ρ}^{(σ + 1)} - D_{z,, μ, ν, ρ}^{(σ)}}{c Δ t} = \frac{a_{- 2} H_{y, μ - 2, ν, ρ}^{(σ + 1 / 2)} + a_{- 1} H_{y, μ - 1, ν, ρ}^{(σ + 1 / 2)} + a_{0} H_{y,, μ, ν, ρ}^{(σ + 1 / 2)} + a_{1} H_{y, μ + 1, ν, ρ}^{(σ + 1 / 2)}}{Δ x} - \\ - \frac{a_{- 2} H_{x, μ, ν - 2, ρ}^{(σ + 1 / 2)} + a_{- 1} H_{x, μ, ν - 1, ρ}^{(σ + 1 / 2)} + a_{0} H_{x,, μ, ν, ρ}^{(σ + 1 / 2)} + a_{1} H_{x, μ, ν + 1, ρ}^{(σ + 1 / 2)}}{Δ y}, \end{array}

In (3а) ‒ (8а):

a_{- 2} = 1 / 6

,

a_{- 1} = - 1

,

a_{0} = 1 / 2

,

a_{1} = 1 / 3

. For the time derivatives we used symmetric (central) approximation.

In order to analyze the errors connected with the numerical dispersion, numerical dissipation and anisotropy of the used difference scheme, let us consider the solution of Maxwell equations on a grid as a plane monochromatic wave:

E_{x, y, z} (μ Δ x, ν Δ y, ρ Δ z, σ Δ t) = E_{x, y, z}^{(0)} \exp [i (k_{x} μ Δ x + k_{y} ν Δ y + k_{z} ρ Δ z - ω σ Δ t)],

H_{x, y, z} (μ Δ x, ν Δ y, ρ Δ z, σ Δ t) = H_{x, y, z}^{(0)} \exp [i (k_{x} μ Δ x + k_{y} ν Δ y + k_{z} ρ Δ z - ω σ Δ t)] .

Here

ω

is a frequency,

k_{x, y, z}

are the Cartesian components of the wave vector

k

. Substituting (16), (17) in (10) ‒ (15), we obtain a system of six linear equations:

(\begin{matrix} R & 0 & 0 & 0 & - Q_{z} & Q_{y} \\ 0 & R & 0 & Q_{z} & 0 & - Q_{x} \\ 0 & 0 & R & - Q_{y} & Q_{x} & 0 \\ 0 & G_{z} & - G_{y} & R & 0 & 0 \\ - G_{z} & 0 & G_{x} & 0 & R & 0 \\ G_{y} & - G_{x} & 0 & 0 & 0 & R \end{matrix}) \cdot (\begin{matrix} E_{x}^{(0)} \\ E_{y}^{(0)} \\ E_{z}^{(0)} \\ H_{x}^{(0)} \\ H_{y}^{(0)} \\ H_{z}^{(0)} \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}),

where

Q_{β} = \exp (- i k_{β} β) \nabla_{β} \exp (i k_{β} β) = [\exp (- 2 i k_{β} Δ β) - 6 \exp (- i k_{β} Δ β) + 3 + 2 \exp (i k_{β} Δ β)] / (6 Δ β)

,

G_{_{β}} = \exp (- i k_{_{β}} β) Δ_{_{β}} \exp (i k_{_{β}} β) = - [\exp (2 i k_{β} Δ β) - 6 \exp (i k_{β} Δ β) + 3 + 2 \exp (- i k_{β} Δ β)] / (6 Δ β)

, and

R = c^{- 1} \exp (- i ω t) \nabla_{t} \exp (i ω t) = (2 i / c Δ t) \sin (ω Δ t / 2)

is symmetric (central) differentiation operator. We can find the dispersion relation considering the condition of existence of solution of system (18) (the equality of determinant of the matrix coefficient in (18) before

E_{x, y, z}^{(0)}

and

H_{x, y, z}^{(0)}

to zero):

\sin^{2} (ω Δ t / 2) / {(c Δ t / 2)}^{2} = F (k_{x}, Δ x) + F (k_{y}, Δ y) + F (k_{z}, Δ z) .

Here, in (12а)

F (k_{β}, Δ β) = [25 + 2 \cos (3 k_{β} Δ β) - 18 \cos (k_{β} Δ β) - 9 \cos (2 k_{β} Δ β)] / [18 {(Δ β)}^{2}] .

Decomposing (19) as a series we obtain the dispersion relation connecting

ω

and

k

:

ω = c k [1 + {(c k Δ t)}^{2} / 24 - (Δ x^{4} k_{x}^{6} + Δ y^{4} k_{y}^{6} + Δ z^{4} k_{z}^{6}) / 30 k^{2} + o (Δ t^{2} + Δ x^{4} + Δ y^{4} + Δ z^{4})] .

From (19) ‒ (21) it follows, that the numerical scheme we use does not possess numerical dissipation (

ω

is a real function of

k

). The second and the third terms in the right part of (21) arise due to the time and space derivatives approximation errors in (1), (2). Their presence leads to the effect of grid dispersion (the dependence of the phase velocity of the monochromatic plane-wave component of the field in the Fourier space (16), (17) on the frequency of this component). Let us especially remark that the third term in right part of (21) not only on the modulus, but also on the direction of the wave-vector, which leads to the anisotropy of a grid. Different signs before the second and the third terms indicate that the errors arising due to the approximation of time and spatial derivatives partially compensate one another.

Let us analyze the stability of our numerical scheme using Von Neumann approach [13]. Let us represent (10) ‒ (15) as

(\begin{matrix} θ_{1} & θ_{2} & θ_{3} & 0 & θ_{4} & θ_{5} \\ θ_{6} & θ_{7} & θ_{8} & - θ_{4} & 0 & θ_{9} \\ θ_{10} & θ_{11} & θ_{12} & - θ_{5} & - θ_{9} & 0 \\ 0 & θ_{13} & θ_{14} & 1 & 0 & 0 \\ - θ_{13} & 0 & θ_{15} & 0 & 1 & 0 \\ - θ_{14} & - θ_{15} & 0 & 0 & 0 & 1 \end{matrix}) \cdot (\begin{matrix} E_{x, μ, ν, ρ}^{(n)} \\ E_{y, μ, ν, ρ}^{(n)} \\ E_{z, μ, ν, ρ}^{(n)} \\ H_{x, μ, ν, ρ}^{(n - 1 / 2)} \\ H_{y, μ, ν, ρ}^{(n - 1 / 2)} \\ H_{z, μ, ν, ρ}^{(n - 1 / 2)} \end{matrix}) = (\begin{matrix} E_{x, μ, ν, ρ}^{(n + 1)} \\ E_{y, μ, ν, ρ}^{(n + 1)} \\ E_{z, μ, ν, ρ}^{(n + 1)} \\ H_{x, μ, ν, ρ}^{(n + 1 / 2)} \\ H_{y, μ, ν, ρ}^{(n + 1 / 2)} \\ H_{z, μ, ν, ρ}^{(n + 1 / 2)} \end{matrix}) .

In (22) the 15 independent elements of matrix

\hat{T}

describing the transfer to the next step are given as follows:

θ_{1} = 1 + (с^{2} Δ t^{2} / ε) (G_{y} Q_{y} + G_{z} Q_{z})

,

θ_{2} = - (с^{2} Δ t^{2} / ε) G_{y} Q_{x}

,

θ_{3} = - (с^{2} Δ t^{2} / ε) G_{z} Q_{x}

,

θ_{4} = - (с Δ t / ε) G_{z}

,

θ_{5} = (с Δ t / ε) G_{y}

,

θ_{6} = - (с^{2} Δ t^{2} / ε) G_{x} Q_{y}

,

θ_{7} = 1 + (с^{2} Δ t^{2} / ε) (G_{x} Q_{x} + G_{z} Q_{z})

,

θ_{8} = - (с^{2} Δ t^{2} / ε) G_{z} Q_{y}

,

θ_{9} = - (с Δ t / ε) G_{x}

,

θ_{10} = - (с^{2} Δ t^{2} / ε) G_{x} Q_{z}

,

θ_{11} = - (с^{2} Δ t^{2} / ε) G_{y} Q_{z}

,

θ_{12} = 1 + (с^{2} Δ t^{2} / ε) (G_{x} Q_{x} + G_{y} Q_{y})

,

θ_{13} = с Δ t Q_{z}

,

θ_{14} = - с Δ t Q_{y}

,

θ_{15} = с Δ t Q_{x}

.

The necessary condition of the scheme stability is the existence of roots $Λ$ of the following characteristic polynomial:

{(Λ - 1)}^{2} {(Λ^{2} - Λ {(c^{2} Δ t^{2} / ε) [G_{x} Q_{x} + G_{y} Q_{y} + G_{z} Q_{z}] + 2} + 1)}^{2} = 0,

which follows from the condition

\det (\hat{T} - Λ \hat{I}) = 0

inside a circle of unit radius on a complex plane (

\hat{I}

is a unit matrix). Using the explicit view of the difference operators

Δ_{β}

and

\nabla_{β}

(back and forward differentiating) we obtain the condition of the scheme stability as a following:

(c^{2} Δ t^{2} / ε) [F (k_{x}, Δ x) + F (k_{y}, Δ y) + F (k_{z}, Δ z)] \leq 4.

The last inequality must be satisfied for any value of the wave vector. Since

\max {(25 + 2 \cos 3 x - 18 \cos x - 9 \cos 2 x) / 18} = 9 / 4

, when we can finally find from (24) the condition on the maximum value of the time step, such, that the numerical scheme remains stable:

Δ t \leq (4 ε^{1 / 2} / 3 c) / {[{(Δ x)}^{- 2} + {(Δ y)}^{- 2} + {(Δ z)}^{- 2}]}^{1 / 2}

.

Acknowledgments

We acknowledge financial support from the Russian Foundation for Basic Research (Grant No. 13-02-00324) and the support from a grant of the President of the Russian Federation for state support of leading scientific schools (Grant No. NSh-3796.2014.2).

References and links

1. J. K. Gansel, M. Thiel, M. S. Rill, M. Decker, K. Bade, V. Saile, G. von Freymann, S. Linden, and M. Wegener, “Gold helix photonic metamaterial as broadband circular polarizer,” Science 325(5947), 1513–1515 (2009). [CrossRef] [PubMed]

2. T. Yoshioka, T. Ogata, T. Nonaka, M. Moritsugu, S.-N. Kim, and S. Kurihara, “Reversible-photon-mode full-color display by means of photochemical modulation of a helically cholesteric structure,” Adv. Mater. 17(10), 1226–1229 (2005). [CrossRef]

3. G. De Filpo, F. P. Nicoletta, and G. Chidichimo, “Cholesteric Emulsions for Colored Displays,” Adv. Mater. 17(9), 1150–1152 (2005). [CrossRef]

4. K. Claborn, E. Puklin-Faucher, M. Kurimoto, W. Kaminsky, and B. Kahr, “Circular dichroism imaging microscopy: application to enantiomorphous twinning in biaxial crystals of 1,8-dihydroxyanthraquinone,” J. Am. Chem. Soc. 125(48), 14825–14831 (2003). [CrossRef] [PubMed]

5. Z. Yang, M. Zhao, and P. Lu, “How to improve the signal-to-noise ratio for circular polarizers consisting of helical metamaterials?” Opt. Express 19(5), 4255–4260 (2011). [CrossRef] [PubMed]

6. Z. Y. Yang, M. Zhao, P. X. Lu, and Y. F. Lu, “Ultrabroadband optical circular polarizers consisting of double-helical nanowire structures,” Opt. Lett. 35(15), 2588–2590 (2010). [CrossRef] [PubMed]

7. M. Thiel, M. Decker, M. Deubel, M. Wegener, S. Linden, and G. von Freymann, “Polarization stop bands in chiral polymeric three-dimensional photonic crystals,” Adv. Mater. 19(2), 207–210 (2007). [CrossRef]

8. N. N. Potravkin, I. A. Perezhogin, and V. A. Makarov, “Numerical solution of Maxwell equations by a finite-difference time-domain method in a medium with frequency and spatial dispersion,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 86(5 Pt 2), 056706 (2012). [CrossRef] [PubMed]

9. G. A. Gryaznov, V. A. Makarov, I. A. Perezhogin, and N. N. Potravkin, “Modeling of nonlinear optical activity in propagation of ultrashort elliptically polarized laser pulses,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 89(1), 013306 (2014). [CrossRef] [PubMed]

10. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Pergamon, 1968).

11. Y. Liu, “Fourier Analysis of Numerical Algorithms for the Maxwell Equations,” J. Comput. Phys. 124(2), 396–416 (1996). [CrossRef]

12. W. H. Press, Numerical Recipes in Fortran: The Art of Scientific Computing, 2nd ed. (Cambridge University, 1992).

13. R. Richtmyer and K. Morton, Difference Methods for Initial-Value Problems (Wiley, 1967).

Ultrashort elliptically polarized laser pulse interaction with helical photonic metamaterial

Abstract

1. Introduction

2. Formulation of the problem

3. The discussion of results

4. Conclusion

5. Appendix

Acknowledgments

References and links

Cited By

Figures (6)

Equations (24)

Optical Materials Express