Wave propagation in metamaterials mimicking the topology of a cosmic string

We study the interference and diffraction of light when it propagates through a metamaterial medium mimicking the spacetime of a cosmic string, -- a topological defect with curvature singularity. The phenomenon may look like a gravitational analogue of the Aharonov-Bohm effect, since the light propagates in a region where the Riemann tensor vanishes being nonetheless affected by the non-zero curvature confined to the string core. We carry out the full-wave numerical simulation of the metamaterial medium and give the analytical interpretation of the results by use of the asymptotic theory of diffraction, which turns out to be in excellent agreement. In particular we show that the main features of wave propagation in a medium with conical singularity can be explained by four-wave interference involving two geometrical-optics and two diffracted waves.


Introduction
Transformation optics [1,2] has become a subject of considerable interest in the recent years given that it provides the ability to control electromagnetic waves in an unprecedented way. The analogy between media and geometry, combined with the advances in the design of structured metamaterials [3], has motivated extensive interest for researchers to develop novel optical devices with unusual properties, some of them inspired by cosmological models, such as omnidirectional light concentrators [4][5][6][7][8][9], optical black holes [10][11][12][13], Hawking radiation in the laboratory [14], optical wormholes [15], optical cosmic strings [16,17], among others. On the other hand, topological defects are of special interest for photonics. They can strongly influence the lightmatter interaction, may induce singularities in the light fields, generate optical beams with orbital angular momenta, produce matter vortices, etc. (see [18,19] and references therein).
In this paper, we address the problem of wave propagation in a metamaterial mimicking the geometry of a topological defect with conical curvature. On the one hand, our study is largely motivated by the corresponding cosmological analogue,the cosmic string, -a topological defect that may have been formed during a phase transition in the early universe [20][21][22]. On the other hand, such defects are known to appear in elastic solids and nematic liquid crystals (called disclinations or wedge dislocations) and they are actively studied [23][24][25][26][27][28]. Recently, it was also pointed out that such media with conical singularities may find important applications in photonics as omnidirectional beam steering elements [16]. The study of this kind of defects in the geometrical-optics limit, or equivalently null geodesics in general relativity, has been performed by different approaches given that analytical expressions can be easily obtained [16,26,29,30]. However, the wave picture is more complex and not many detailed studies are available [31][32][33][34]. Here, we give the full interference and diffraction description of wave propagation in a medium with conical curvature, by providing the full-wave numerical simulation and comparing the results with the analytical theory developed recently [33,34]. Furthermore, we examine different asymptotic approximations -the uniform asymptotic theory [35,36] and the geometrical theory of diffraction [37], -applied for this kind of media with conical singularities. Finally, we study the effects of the wavefront curvature on the results by contrasting the case of the finite-distance source with that for the plane-wave incidence.

Medium parameters
The development of transformation optics in recent years has established the relationship between geometry and material properties [1,2]. In particular, it has been shown that light propagation in a curved vacuum manifold is formally equivalent to light propagation in a medium filled with an inhomogeneous anisotropic material embedded in a flat Minkowski spacetime [38].
To mark out the framework of the model in which we will work, let us consider the spacetime metric in arbitrary coordinates where g ij is the spatial part of the metric tensor with i, j = 1, 2, 3. Given that our goal is to study wave effects in a spatially anisotropic material, for simplicity, we consider the case when space and time components of the metric are decoupled, g 0i = 0. (For the conical geometry we will consider later on, this assumption corresponds to the nonspinning cosmic string.) In this case, there is no coupling between the electric and magnetic fields for the equivalent medium [38] and the constitutive relations take the usual form for an anisotropic material with the permittivity and permeability tensors given in Cartesian coordinates by [39] Here, g ij is the inverse of g ij , and g is the determinant of the full spacetime metric. In this way, the spacetime geometry is mapped into the parameters of a medium. Next, consider the conical geometry which can be described in cylindrical coordinates by the line element where β is a conical parameter which corresponds to the removal (0 < β < 1) or insertion (β > 1) of a wedge with an angle 2∆. We define ∆ ≡ π(1 − β), given that the removal of the wedge (β < 1) will be the only case studied in this work, therefore ∆ > 0 is assumed. The metric (4) with β = 1 − 4µ S describes the spacetime of a static gravitating cosmic string [22] with a mass per unit length µ S confined at the origin r = 0 along the zaxis (by using units in which G = c = 1). It also represents the effective geometry of a linear topological defect in condensed-matter systems (e.g., disclination in nematic liquid crystals [25,26], or wedge dislocation in elastic solids [23]). It should be noted that cosmic strings formed in the early universe are expected to have a rather small angular deficit, typically with ∆ ∼ 10 −7 [21,22]. Here, we provide simulation results for the larger scale range, 0 < ∆ π/2, which is of broader interest in applications and can be easily visualized. Yet, our analytical results are valid in the limit ∆ 1 as well. By applying the transformation (3) to the conical geometry (4), one can find the medium parameter tensors to be diagonal in cylindrical coordinates (r, ϕ, z), which can also be written in the usual Cartesian base (x, y, z) as Note that the parameters (6) vary with the angle ϕ but are independent of r, z. This fact is clearly seen in figure 1 where the spatial distributions of the tensor are plotted. Another distinctive feature is the singularity at the origin, which we will discuss in detail later on. For the case of a spinning cosmic string, the medium parameters can be found in Ref. [40].

Curvature singularity
A conical spacetime (4) is an example of a geometry with a singularity at which the curvature cannot be calculated by standard formulas [41,42]. Yet one can determine the integral of the curvature by applying the Gauss-Bonnet theorem that relates the Gaussian curvature K over some region Σ to the geodesic curvature k g calculated along the smoothly closed boundary ∂Σ [43]: On the other hand, the rhs of equation (7) is the angle χ a transported frame is rotated through as a result of parallel propagation around ∂Σ [42]: withχ being the total geodesic curvature. Intuitively, the geodesic curvature measures how far a curve is from being a geodesic, the shortest path between two points [43].
To calculate the curvature for the metric (4), let us determineχ over a curve which encloses the singularity. We take Σ to be a disk of radius R, at z = 0 plane, centred at the singularity r = 0, and its circumference ∂Σ to be a curve γ parametrized by arc-length s [see figure 2(a)]. We obtaiñ with Γ r ϕϕ being the relevant Christoffel symbol andφ ≡ dϕ/ds is the angular velocity in units of the arc-length. For the metric (4), Γ r ϕϕ = −β 2 R andφ = −1/(βR), so that the product Γ r ϕϕφ = β is constant. That givesχ = 2πβ and the frame rotation which is also equal to the angular deficit 2∆. Hence, the integral curvature of a disk centred at the singularity is 2π(1 − β) independently of its radius R. Next, we consider the closed path that does not enclose the singularity, as shown in figure 2(b). Since the product Γ r ϕϕφ = β calculated on the arc is independent of the radius, the integrals over γ 1 and γ 3 , being in the opposite directions, compensate each other. The paths γ 2 and γ 4 are geodesics, therefore they do not contribute to the geodesic curvature. As a result, for the total path we getχ = 0. To apply the Gauss-Bonnet theorem for this case, we have to take into account that, instead of a smoothly closed curve, we have a piecewise curvilinear polygon. The theorem still holds, but we must include a correction for the vertices of the polygon where the curve is not smooth [43]. We obtain χ = 2π − 4 · (π/2) = 0.
Summarizing the above arguments for the two loops, one can conclude that the Gaussian curvature K for the conical spacetime (4) can be defined in terms of the δ-function with its singularity at the origin. Moreover, in a 2-dimensional manifold, only one component of the Riemannian tensor is independent and can be related to the Gaussian curvature as K = R rϕ rϕ [43]. Consequently, the corresponding components of the Riemannian and Ricci tensors can be found in the form [41,42] where δ 2 (r) is a two-dimensional Dirac δ-function in flat space, which is equal to δ(x)δ(y) in Cartesian coordinates. In the following, we only consider the case β < 1 that corresponds to a singularity with positive curvature. The value of β = 1 gives an everywhere flat Minkowskian spacetime with no curvature. It should be noted, that the conical singularity is topological in nature and cannot be eliminated by a coordinate transformation. Therefore, when we consider a metamaterial analogue in flat spacetime by applying the transformation (3) to the conical geometry (4), the singularity should persist, being embedded, in this case, in the medium parameters. We also notice an interesting property which follows from the Gauss-Bonnet theorem: for the conical spacetime no curve enclosing the singularity could be a geodesic line since it must have a zero geodesic curvature which is in contradiction with equation (10).

Wave equation
Having found the medium parameters which mimic the conical topology of a cosmic string, we now consider the propagation of a monochromatic electromagnetic wave with frequency ω in such a medium. The Maxwell equations supplemented by the constitutive relations (2) lead to the wave equation for a time-harmonic electric field vector E ω in an anisotropic medium [44] where ε and µ −1 denote the permittivity and inverse permeability tensors, respectively. Similar to Ref. [11], we consider a transverse electric (TE) wave propagating at z = 0 plane, E ω (x, y) = E(x, y)ẑ, where E is the out-of-plane z component of the electric field. Then, the corresponding wave equation in the two-dimensional (x, y) space is given by where µ xx , µ xy , µ yy are the components of the inverse permeability tensor µ −1 and ε zz corresponds to the permittivity ε. A simple inspection of Eq. (13) shows that the medium parameters can further be simplified. Indeed, since only one component of ε has entered into the equation, and it is constant [see Eq. (6)], the permittivity can be chosen to be isotropic with all its diagonal elements equal to ε zz and the off-diagonal terms equal to zero. In this case, only the magnetic anisotropy is left. Moreover, one can redefineẼ = βE, to obtaiñ where δ ij = diag (1, 1, 1). In this way, the metamaterial has the permittivityε of free space and the permeabilityμ inhomogeneous in two dimensions, simplifying its design. Finally, if we write the wave equation for this metamaterial in polar coordinates, we obtain which is precisely the equation used in Refs. [31][32][33][34] for the cosmic string background (4).
In the next section, we will solve numerically the wave equation (13) for a TEpolarized wave propagating in a magnetically anisotropic medium. It should be noted, that the transverse magnetic (TM) wave can also be considered in a similar medium after the substitutions: ε ij µ ij , E → H. In such a case, the medium is electrically anisotropic, and one should solve the equation for the magnetic field H(x, y).

Numerical simulation
To solve numerically Eq. (13) we use a rectangular 2D geometry of (x, y) space as shown in figure 3. We consider a 2D wave produced by a point source on the left boundary at a distance r 0 from the string S. The source is an electric line current perpendicular to the domain which generates TE polarized cylindrical waves. The wave equation is solved for the electric field E. The distances are normalized by the wavelength λ, which makes the equation independent of the frequency. Therefore, it can be solved for any wavelength as long as it is larger than the unit cells of the metamaterial and within its operating bandwidth [2]. The same value for the source-string distance r 0 = 10λ is used for all the simulations in this work. The simulation was carried out by using the COMSOL Multiphysics software package and the figures were processed in Python.
In figure 4 we show typical results of our simulation for the electric field E in a medium which mimicks cosmic-string topology. One can see that the emitted wave propagating around the string interferes with itself giving rise to a characteristic pattern. This self-interference is caused by the conical curvature (singularity) at the string location. By plotting the squared electric field norm in Figure 4(b), we also see that the wave field is amplified in some areas behind the string due to diffraction. These phenomena will be discussed in detail in the next section along with the quantitative analysis of the results. In order to quantify the effect that the string causes on the wave field, we will also introduce the amplification factor F ≡ |E/E 0 |, where the electric field E is normalized by its value E 0 in the absence of the string.

Geometrical optics limit
First of all, our goal is to understand how the interference pattern observed in numerics is formed. We begin with a simple situation: consider a point source A emitting light beams in the neighbourhood of a string S as shown in figure 5(a). The conical singularity has the effect of focusing geodesics whenever β < 1. This can be readily seen by changing the angular coordinate to θ = βϕ, where θ spans the range 2πβ = 2π − 2∆ in the Minkowski space with a wedge removed and the faces of the wedge identified [22]. To perform the angular transformation, we choose to place the cut edge along the line AS linking the source and the string. In this way, the source A is doubled to A − and A + , each on a different side of the wedge, as depicted in figure 5(b) (see [34] for further details). In the new representation, the geodesics are simply straight lines and it is easy to see that the region −∆ < θ < ∆, called the double-imaging region, is illuminated by both image sources A − and A + . Indeed, an observer located at P 1 in figure 5 will see two images of the source A. On the other hand, an observer at P 2 will only see one image, being beyond the shadow line for the second one. Thus, the conical singularity caused by the string gives rise to a double-imaging region behind it, in which two image sources interfere. From the view of figure 5(a), there is only one source but the beams, passing on opposite sides of the string, intersect producing an interference pattern along and across the line of sight.
The beam interference can be described in the most simple terms as the geometrical optics (GO) solution. In this limit, the field will be the sum of two contributions [34] where k is the wavenumber in free space, s ± = r 2 + r 2 0 + 2rr 0 cos(∆ ± βϕ) are the path lengths of the GO waves from the image sources A + and A − to an observation point [in (r, θ) space they are cylindrical waves for our geometry], and h ± ≡ H(∆ ± βϕ), with H(x) being the Heaviside step function that guarantees that the double-imaging range is determined by −∆/β < ϕ < ∆/β. Equation (16) allows to understand the appearance of bright and dark lines behind the string (which we will refer to as "antinodal" and "nodal" lines, respectively) seen in figure 4. Indeed, to have constructive or destructive interference, disregarding the slowly varying pre exponential factors, the path length difference, s − − s + , must be equal to an integer number of half wavelengths In figure 6 these lines are superimposed on the wave pattern (shown in color) simulated by solving the wave equation numerically. A good agreement is observed within the double-imaging region. However, one cannot explain the inerference pattern outside that region and the modulation of the intensity over the antinodal lines within the GO approximation.

Geometrical theory of diffraction
To advance further, we use the geometrical theory of diffraction (GTD) [37] by considering an additional contribution caused by diffraction. We write the wave field E at a point (r, ϕ) as the sum of the two GO waves given by equation (16) and two diffracted (D) waves where D ± are the diffraction coefficients (or directivity functions) given by [34] The physical interpretaion of the D waves is the following: they are the waves that go from the source to an observer but hitting the string, following the shortest path [33,34].
The coefficients D ± give the amplitude of the D wave as a function of the direction ϕ and they are singular at ϕ = ±∆/β, the boundary of the double-imaging region. In figure 7 we compare the results of numerical simulation in the metamaterial with the GTD solution (18). It is seen that the wave effects behind the string are very well reproduced except at the boundary of the double-imaging region, where the Dwave terms diverge. We also observe that the double-imaging effect results in the field amplification behind the string and, due to diffraction, the factor F can even be greater than 2.
It should be emphasized that the D-wave terms are not the next order terms in the expansion over the wavelength λ, as it would be in the case of a diffraction  Figure 7. Comparison of the wave patterns obtained from (a) metamaterial simulation, (b) geometrical theory of diffraction (GTD); for a medium with ∆ = π/5. on a compact object. For the conical defect, the relevant parameter is the Fresnel number N F ∼ r∆ 2 /λ, where r is the distance from the string to the observer and ∆ is proportional to the mass of the string (the singularity strength). For ∆ and λ fixed, one can always find the distance r, for which N F ∼ 1, and the D and GO waves are therefore of the same order of magnitude [34]. The GO terms dominate whenever N F 1. The four-wave interference description of the GTD may also be helpful to explain the intensity modulation over the antinodal lines. Both D waves follow the path of r+r 0 , therefore the path difference between the GO and D waves will be r + r 0 − s ± . Also, we have to take into account that the D wave acquires a phase shift of 3π/4 by hitting the string, as seen in the diffraction coefficients (19). Therefore, we would expect the maxima and minima to occur when these two conditions are fulfilled simultaneously: with m, n being non-negative integers: 0, 1, 2, . . .. These constant-phase lines are hyperbolas in (r, θ) space [34]. It is easy to see that, by subtracting equations (20) and (21), one obtains equation (17) with q = n − m. That means that the intersection points of the hyperbolas lie on the nodal and antinodal lines (see figure 8). The maxima are found whenever both indices n and m are even, while the saddle points (local minima on the antinodal lines) correspond to simultaneously odd n and m. The interference points can be labelled with the indices (n, m) and be associated to "Fresnel observation zones" (see [34] for further details). Summarizing, the main features of the wave pattern observed in numerics can be explained as the interference of four characteristic waves of the GTD. It is known from singular optics [45,46] that when three or more waves interfere in two dimensions, the intensity vanishes at points, rather than on lines, that may produce optical singularities. This is not the case for the conical space and the type of waves we consider, where four characteristic waves (two GO and two D waves) have only two independent phase differences, due to the symmetry, hence the intensity vanishes on nodal lines.

Uniform asymptotic theory
The GTD we have used to explain the numerical results is quite satisfactory in the whole region of interest except at the boundary of the double-imaging region where the solution diverges. This nonuniform solution can however be "regularized" eliminating the singularity in the framework of the uniform asymptotic theory (UAT) [35,36]. We write the solution for the wave field in the form [34] which is the sum of the penumbra field defined in terms of the Fresnel integral F(w) = e −iπ/4 π −1/2´w −∞ e iu 2 du and the diffracted field with modified diffraction coefficients Here, we use the notations w ± = σ ± k(r + r 0 − s ± ) and σ ± = sgn(∆ ± βϕ). Note that the solution (22) is uniform, since it is finite and continuous across the boundary of the double-imaging region. It can be shown [34], that far from the boundary (w ± 1), both asymptotics, the GTD (18) and UAT (22), coincide. In figure 9 we compare the  Figure 9. Comparison of the wave patterns obtained from: (a) metamaterial simulation, (b) uniform asymptotic theory (UAT); for a medium with ∆ = π/5. results obtained by numerical simulation in a medium with the UAT solution given by equation (22). We see an excellent agreement in the entire region.
The numerical results for another value of ∆ = π/4 are shown in figure 10 in comparison with the analytical theories discussed above. This figure allows to compare all types of approximations in one view and observe the corrections that each one gives. With the GO approximation [ figure 10(b)] one can already observe the main interference effects in the double-imaging region. Nevertheless, there is an abrupt change of the field at the boundary between the single-and double-imaging zones with no wave effects in the former. The wave pattern is substantially improved by use of the diffraction terms in the GTD, which also add the modulation in the antinodal lines [figure 10(c)]. Lastly, the UAT [figure 10(d)] provides a smooth transition between the single-and double-imaging regions, having the least discrepancy with the numerics. How small the discrepancy is, one can appreciate from 1D plots of figure 11, where the amplification F is shown as a function of ϕ for two fixed distances r. One can see that the agreement

Plane wave propagation in conical space
If the distance from the source to the string goes to infinity, r 0 → ∞, one can neglect the effects of the wavefront curvature and consider the incidence of a plane wave. In this limit, the field at a point (r, ϕ) is given by [33] E = e ikr cos(∆+βϕ) F(u + ) + e ikr cos(∆−βϕ) F(u − ) (24) with u ± = √ 2kr sin[(∆ ± βϕ)/2]. We did not carry out any numerical simulations for this limiting case, however, it would be interesting to compare the analytical results (24) with the similar results for a finite-distant source. For a plane wave, the nodal and antinodal lines should be straight lines in (r, θ) space, parallel to the line of sight with q being an integer. After the angular transformation θ = βϕ, those lines will still be almost parallel to the x-axis, as long as sin(βϕ) ≈ βϕ for small arguments. This is seen in figure 12(a) where the spatial distribution of the amplification factor F (equal to |E| 2 in this case) is presented for ∆ = π/6. The typical spacing between the fringes along the y-axis, as follows from equation (25), is approximately δy ≈ λ/(2β sin ∆), and it is independent of the distance from the string. For plane waves, the hyperbolic lines of constant phase between the GO and D waves become parabolic, r[1−cos(∆±θ)] = const, and the intensity maxima may be found at the intersection points of the curves with m, n being non-negative integers, as it is depicted in figure 12(b).
Given an excellent correspondence with the numerical results, one can use equations (26) to predict the location of some typical maxima. For instance, the position of the highest diffraction maximum can be estimated by substituting m = n = 0 and ϕ = 0. For ∆ = π/6, this gives x ≈ 2.8λ. The next-order maxima (n = 2 and m = 0 and vice versa) are at x ≈ 5.9λ and y ≈ ±1.2λ, while the fringe separation for this ∆ should be δy ≈ 1.2λ. All these values are seen to be in agreement with figure 12.

Concluding remarks
Transformation optics is a powerful tool which allows one to design an artificial metamaterial medium in which the light behaves in a similar way as if it were in the curved space. In this paper, we have analyzed the medium properties and the propagation of electromagnetic waves through the effective media which mimic the conical topology of a static cosmic string -a 1D topological defect. A conical spacetime is an example of a geometry which has a singular δ-like curvature at the string core. This singularity has important physical consequences similar to those which arise in the Aharonov-Bohm setting. Photons or particles constrained to move in a region where the Riemann tensor vanishes may nonetheless exhibit physical effects arising from non-zero curvature confined to the string core, i.e., a region from which they are excluded. This is a gravitational analogue of the Aharonov-Bohm effect for a static mass distribution [47][48][49]. The physical manifestations of this effect can be summarized as follows: (i) In the geometrical-optics limit, the light beams or photons propagating on different sides of the string can intersect, and their relative deflection cannot, in general, be transformed away. This is analogous to the electromagnetic Aharonov-Bohm effect, where the vector potential can be transformed to zero only in regions containing no closed loops enclosing the solenoid [47][48][49].
(ii) The deflection angle ∼ ∆ is independent of the impact parameter, but depends on the strength of the singularity, i.e., on the mass of the string [21].
(iii) The beams intersection gives rise to a double imaging effect.
(iv) In the double-imaging region, the interference pattern should appear with the spacing between the fringes asymptotically ∼ λ/(2β sin ∆), which is independent of the distance.
(v) The diffraction effects are manifested in the intensity modulation along the antinodal lines and in the wave pattern beyond the geometrical-optics domain.
(vi) The Fresnel number N F ∼ r∆ 2 /λ depends on the radial coordinate, hence, the Fresnel observation zones can be introduced associated with the intensity amplification caused by diffraction [34].
In the present paper, we are mainly concerned with the wave effects (iv) and (v). We carry out the full-wave numerical simulation of the metamaterial medium which mimics the cosmic-string topology. The obtained results are interpreted by use of the asymptotic theories of diffraction which, thus far, have been essentially applied to obstacles with clearly defined boundaries (a half-plane, a slit, etc. [36]). Here, we apply them to a topological defect and obtain the excellent agreement with the numerical results. We notice the advantage of the asymptotic theories. The existing solutions for similar problems in the form of integral representation [31] or infinite series [32] exclude rigorous analysis necessary for practical applications. The infinite series solution, for instance, is poorly convergent for the distances greater than about a wavelength. These limitations are overcome by the GTD and UAT we have used in this work, since they capture the main physics and a few terms are only needed to achieve the required accuracy. At first sight, the emergence of two image sources in conical space followed by their mutual interference, may look like a familiar Young's double-slit interference experiment. Our results, however, indicate that the two models are conceptually different. The Young's double-slit setting has a characteristic length -the distance between the slits d, that implies that the interference fringe scales as ∼ λ r/d [50]. On the contrary, for the conical space, there is no characteristic length but the deficit angle ∆. As a result, the interference and diffraction effects manifest differently from those in the Young's experiment (see items (iv)-(vi) above). All of these peculiar features open up new opportunities for applications in photonics and other fields. For instance, such metamaterial media with conical singularities can act as omnidirectional beam steering devices [16], beam splitters [16,51,52], diffraction-control elements [17], etc. Our results can also potentially be applied to light [26,53] or sound propagation [25,28] near a linear topological defect in nematic liquid crystals. It would also be interesting to extend our approach to other geometries or types of topological defects.