Spontaneous isotropy breaking for vortices in nonlinear left-handed metamaterials

TL;DR

This paper investigates how vortex beams in nonlinear left-handed metamaterials spontaneously break isotropy, forming polygonal patterns with symmetries depending on the vortex charge, explained through an effective Landau-Ginzburg model.

Contribution

It introduces a novel theoretical framework linking symmetry breaking in vortex beams to nonlinear metamaterials using an effective field theory approach.

Findings

Vortex beams induce polygonal symmetry patterns in metamaterials.

The symmetry depends on the vortex charge and material response.

Loop corrections improve the accuracy of theoretical predictions.

Abstract

We explore numerically and analytically the pattern formation and symmetry breaking of beams propagating through left-handed (negative) nonlinear metamaterials. When the input beam is a vortex with topological charge (winding number) $Q$, the initially circular (isotropic) beam acquires the symmetry of a polygon with $Q$, $2Q$ or $3Q$ sides, depending on the details of the response functions of the material. Within an effective field-theory model, this phenomenon turns out to be a case of spontaneous dynamical symmetry breaking described by a Landau-Ginzburg functional. Complex nonlinear dependence of the magnetic permittivity on the magnetic field of the beam plays a central role, as it introduces branch cuts in the mean-field solution, and permutations among different branches give rise to discrete symmetries of the patterns. By considering loop corrections in the effective Landau-Ginzburg field theory we obtain reasonably accurate predictions of the numerical results.

Equations

$$ϵ(E,E^{†})$$

$$\mu (H,H^{†})$$

$$|H{|}^2=\alpha A^2\frac{(1-X^2)\left[{\left(X^2-{\Omega }^2\right)}^2+{\Omega }^2{\gamma }^2\right]}{X^6},$$

$$-\frac{ı}{b}{\partial }_{z}E$$

$$-\frac{ı}{b}{\partial }_{z}H$$

$$\mathcal{L}$$

$${\mathcal{L}}_E$$

$${\mathcal{L}}_H$$

$$-\frac{ı}{b}{\partial }_{z}E$$

$$-\frac{ı}{b}{\partial }_{z}H$$

$${\Phi }_H$$

$${\Phi }_E$$

$${\Phi }_E=\frac{{\mu }^{\prime }}{\mu }{\nabla }_{\perp }\left(\frac{E{\nabla }_{\perp }{E}^{†}-E^{†}{\nabla }_{\perp }E}{\mu }\right).$$

$$E=E_{\mathrm{reg}}+E_{\mathrm{vort}},$$

$$E_{\mathrm{vort}}(r,\varphi ,z)=Z_E(z)R_E(r)e^{ıQ\varphi -ı\Phi (\varphi )},$$

$${\Phi }^{\prime \prime }-ı{\left({\Phi }^{\prime }\right)}^2+2ıQ{\Phi }^{\prime }+ıl^2=0,$$

$$\Phi (\varphi )=\cos \left(\sqrt{Q^2+l^2}\varphi +C_l\right).$$

$$R_E^{\prime \prime }+\frac{1}{r}{R}_E^{\prime }+\left(\frac{\lambda }{r^2}-k^2+{ϵ}_{D0}{\tilde{\omega }}^2\right)R_E+\frac{\alpha {\mu }_c{\tilde{\omega }}^2}{E_c^2}{R}_E^3=0,$$

$$R_E(r)$$

$$Q_l$$

$$G_{E,H}(\mathbf{u})=\frac{2\pi }{\Gamma (Q/2)}\frac{e^{ıQ(\pi /2+\varphi )}}{au}\left(\frac{\sin \left(\left(u-a\right)\Lambda \right)}{u-a}-\frac{\cos \left(\left(u+a\right)\Lambda -\pi Q\right)}{u+a}\right)+$$

$$+\frac{2\pi }{\Gamma (Q/2)}\frac{e^{-ıQ(\pi /2+\varphi )}}{au}\left(\frac{\cos \left(\left(u-a\right)\Lambda +\pi Q\right)}{u-a}-\frac{\cos \left(\left(u+a\right)\Lambda \right)}{u+a}\right).$$

$${\left({\omega }_{0NL}^{(1)}\right)}^2+{\left({\omega }_{0NL}^{(2)}\right)}^2+{\left({\omega }_{0NL}^{(3)}\right)}^2=\frac{1+2{\Omega }^2}{1+|H{|}^2/\alpha E_c^2},$$

$$\delta \mathcal{L}=g_{2,0,0}|{\nabla }_{\perp }E{|}^2+g_{0,2,0}|E{|}^2+g_{0,2,2}|E{|}^2|H{|}^2+g_{2,0,2}|{\nabla }_{\perp }E{|}^2|H{|}^2.$$

$${\Sigma }_E^{(1)}$$

$${\Sigma }_E^{(2)}$$

$$E_{\mathrm{vort}}(r,\varphi ,z)$$

$$\mathcal{L}=\frac{{\left(E_0^{\prime }\right)}^2+\frac{{\left({\Theta }^{\prime }\right)}^2}{r^2}+k^2{E}_0^2}{\mu }+\frac{{\left(H_0^{\prime }\right)}^2+\frac{{\left({\Theta }^{\prime }\right)}^2}{r^2}+k^2{H}_0^2}{ϵ}.$$

$$\frac{ϵ\left({\nabla }_{\perp }^2-H\right)E-{ϵ}^{\prime }{\nabla }_{\perp }E\cdot {\nabla }_{\perp }H}{{ϵ}^2}+\frac{{\mu }^{\prime }}{{\mu }^2}|H{|}^{-1/3}=0,$$

$$\frac{{\left({\Theta }^{\prime }\right)}^2\left(1-\frac{{ϵ}^{\prime }}{ϵ}\frac{E_0}{H_0}\right)-k}{ϵ}+\frac{2{\mu }^{\prime }}{3{\mu }^2}|H{|}^{-1/3}=0.$$

$$ı({\Theta }_0+2\pi n/3)=K_n\log \left[\frac{k\left(1-\frac{{ϵ}^{\prime }}{ϵ}\frac{E_0}{H_0}\right)}{E_0^2+H_0^2}{\sec }^2\left(\frac{Q}{2}\varphi \right)\right],$$

$$\nabla \cdot \hat{D}=\rho =0,\nabla \cdot \hat{B}=0,\nabla \times \hat{E}=-{\partial }_t\hat{B},\nabla \times \hat{H}=4\pi \hat{j}+{\partial }_t\hat{D}={\partial }_t\hat{D}.$$

$$\nabla \times \nabla \times \hat{H}=-{\nabla }^2\hat{H}+\nabla \left(\nabla \cdot \frac{\hat{B}}{\mu }\right)=-{\nabla }^2\hat{H}+\nabla \left(\frac{1}{\mu }\nabla \cdot \hat{B}\right)-\nabla \left(\frac{\nabla \mu }{{\mu }^2}\cdot \hat{B}\right)=$$

$$=-{\nabla }^2\hat{H}+\nabla \left(\frac{1}{\mu }\nabla \cdot \hat{B}\right)-\nabla \cdot \left(\frac{\nabla \mu }{{\mu }^2}\right)\hat{B}-\frac{\nabla \mu }{{\mu }^2}\nabla \cdot \hat{B}=-{\nabla }^2\hat{H}+0+O(1/l^2)+0=-{\nabla }^2\hat{H},$$

$$\nabla \times \nabla \times \hat{H}=\nabla \times ({\partial }_t\hat{D})=-ı\omega \nabla \times \hat{D}=-ı\omega \nabla \times (ϵ\hat{E})=-ı\omega (\nabla ϵ)\hat{E}-ı\omega ϵ\nabla \times \hat{E}=$$

Full text

I Introduction

The idea of a material with negative refraction index was first considered long before it could be realized in experiment, in the now famous paper by Veselago veselago , in 1968. He considered a material with negative electric permeability $\epsilon$ and magnetic permittivity $\mu$, and predicted a number of interesting properties in such systems, among them negative refraction. Only much later did it become possible to combine elements with negative $\epsilon$ and negative $\mu$ at a microscopic level, as a composite metamaterial. First experimental realizations were reported in first ; exper . Negativity, or left-handedness, is typically only achieved in a narrow frequency range, close to the resonant frequency of the conductive elements of the metamaterial. This was the original motivation for studying nonlinear effects in these systems. Nonlinearities can be strengthened by appropriate design at the microscopic level. The study of nonlinear phenomena in metamaterials started with zarov03 . This has become a broad and important field in metamaterials research rmp . Nonlinear phenomena like solitons zarova05 ; sadriv05 , nonlinear surface waves sadriv04 , modulational instability wen ; walasik and ultrashort pulses tsitsas were observed. Other work in left-handed metamaterials relevant for our paper includes, among others, sadriv03 ; kivstamm ; ma ; fan ; korotki ; sarma ; kevrek ; ciral ; symmbreak . We have no intention of being exhaustive in this short review of the literature; we merely mention the results we have directly used or found particularly inspiring.

The focus of our work is the dynamics of symmetry breaking in intensity patterns of electromagnetic waves propagating through a left-handed nonlinear metamaterial. Numerical solutions of the equations of motion reveal that circular (usually Gaussian) input beams turn into polygonal patterns, with some discrete symmetry. This fits the textbook notion of symmetry breaking, more specifically dynamical symmetry breaking. The general theory of dynamical criticality is by now well-developed cross and has been applied to numerous systems rabin . In rabin , a systematic theory of isotropy-breaking transition is presented, though mainly for periodic and quasiperiodic structures (convection in fluids, fluctuations in quasicrystals). The basic mechanism is that the system develops momentum eigenmodes of fixed module but with multiple discrete directions on the sphere $|\mathbf{k}|=\mathrm{const.}$ in momentum space. In nonlinear negative materials, the situation is complicated by the strong frequency dependence of the magnetic permittivity but the same basic logic remains. At a fundamental level, this situation can be understood from the viewpoint of a spatially non-uniform Landau-Ginzburg theory. Quantitative accuracy is however hard to achieve; this requires cumbersome perturbative calculations. Ultimately, numerical work is the best way to describe the patterns in detail; they look like polygons or, occasionally, necklaces, with $C_{3Q}$, $C_{2Q}$ or $C_{Q}$ symmetry, depending on the parameter regime; here, $Q$ is the topological charge of the beam, a property we will discuss in detail in the next paragraph. The paper walasik is very important in this context: it starts from the model derived in wen and studies mainly necklace configurations, which consist of discrete beads (spots of high intensity) distributed more or less uniformly along a circle. The authors find the same $C_{3Q}$ symmetry that we see. Our goal is to gain a detailed understanding of the phenomenon, and move beyond single beams toward collective behavior and interactions.

We have chosen to study this phenomenon on vortices, topologically nontrivial solutions where the phase of the complex electric and magnetic field winds one or more times along a closed line encircling the vortex core. Vortices appear in many systems described by a complex field, i.e., a field with $U(1)$ phase invariance vortbook ; kleinert . In optics, this is just the complex beam envelope of the electric and magnetic field. The phase of any complex wavefunction or field can wind along some closed line around a defect, forming a vortex. Famously, vortices may coexist with the superconducting order ($U(1)$ symmetry breaking) in type-II superconductors or they may exist only in the normal phase, upon destroying the superconductivity (type I). Pattern-forming systems like fluids and soft matter often have rich vortex dynamics rabin . Other examples of vortex matter in nature arise in liquid helium liqhel , Bose-Einstein condensates bec and magnetic systems vortcupsach . In two spatial dimensions, interactions among the vortices lead to a vortex unbinding phase transition of infinite order found by Berezinsky, Kosterlitz and Thouless for the planar XY model kosterthoul . We study a three-dimensional metamaterial but with an elongated geometry, so we treat it as a $2+1$-dimensional system (the $x$ and $y$ coordinates are spatial dimensions and the $z$-direction has the formal role of time). We therefore have a similar situation to the XY model but with different equations of motion and different phenomena.

In addition to direct numerical and analytical study of the equations of motion, we also propose an effective field-theory Lagrangian which gives slightly different equations but captures the key properties of the system. The Lagrangian form makes it easier to understand some of the phenomenology we find in numerical simulations; the foundations of the symmetry breaking are obtained from this model in a natural way. Numerical work is done with original equations of motion, as they are directly grounded in the microscopic physics. The Lagrangian is just a phenomenological tool to facilitate the theoretical understanding. It is difficult (and perhaps impossible) to package the exact original equations in a Lagrangian form because the system is strongly nonlinear and dissipative. Dissipative systems can be encapsulated in a Lagrangian (our Lagrangian is also dissipative) but with some limitations, and there is certainly no general method to write down a Lagrangian for a broad class of dissipative systems.

The structure of the paper is the following. In the next section we describe the model of a nonlinear left-handed metamaterial, following closely the wave propagation equations used in previous research, e.g. in zarov03 ; zarova05 ; sadriv05 and others, which correspond to a specific experimentally realizable metamaterial. We also formulate and motivate the field-theory model of the system. In the third section, we describe our numerical findings, above all the anisotropy of the intensity patterns. The fourth section offers the theoretical explanation for the patterns: first by a direct approximate solution of the propagation equations, and then also from field theory, which makes the physics of the symmetry breaking particularly clear. In the fifth section we briefly discuss how to check our predictions in experiment and how prominent the effects of symmetry breaking are compared to other possible instabilities in realistic metamaterials. The last section sums up the conclusions. We have included some long calculations in the Appendices.

II Wave equations in a nonlinear left-handed metamaterial

We adopt the model of zarov03 ; sadriv05 to describe a left-handed metamaterial with a nonlinear response. Microscopically, the material is realized as a lattice of split-ring resonators and wires. In the terahertz range, this is an experimentally well-studied system exper . In zarova05 , a detailed microscopic derivation is given, starting from the current transport equations in the resonator-wire system. The outcome is a nonlinearity similar to that postulated phenomenologically in zarov03 . We adopt essentially the same model, described by the electric permeability $\epsilon$ and the magnetic permittivity $\mu$:

[TABLE]

with $\alpha=1$ or $\alpha=-1$ for self-focusing or self-defocusing nonlinearity, respectively. Frequency is denoted by $\omega$ and $\epsilon_{D0}$ is the linear part of the permittivity. By $F$, $\gamma$ and $\Gamma$ we denote the filling factor of the material and the electric and magnetic damping coefficients. The equations (1-2) allow us to model also the real (lossless) dielectric response by putting $\gamma=0$. For the magnetic field, the permittivity will in general stay complex even for $\Gamma=0$, as the nonlinear frequency of the resonator rings $\omega_{0NL}$ can always have a nonzero imaginary part. This frequency is related to the magnetic field through the relation ($X\equiv\omega_{0NL}/\omega_{0}$):

[TABLE]

where $\Omega\equiv\omega/\omega_{0}$, $\omega_{0}$ is the eigenfrequency of the rings, and $A$ is a parameter which can be derived microscopically zarova05 ; zarov03 ; sadriv05 ; for our purposes, it can be treated just as a phenomenological parameter. This cubic equation yields three branches for $\omega_{0NL}^{2}$. All these branches are physical and correspond to different possible nonlinear oscillations sadriv05 . Now the equations of motion are just the Maxwell equations in a medium described by (1-2), in the approximation of slowly-changing beam envelopes. We assume an elongated (cylindrical or parallelopipedal) slab of metamaterial, so we can employ the paraxial beam approximation (e.g. book ). The beam is initially collimated along the longitudinal axis $z$ and focuses or defocuses slowly in the transverse $x-y$ plane due to the nonlinearity of the material. The electric and magnetic field $\hat{E}(t;x,y,z),\hat{H}(t;x,y,z)$ are directed along the $z$-axis. From now on, the speed of light is put to unity $c=1$. All the steps in deriving the nonlinear-Schrödinger-like equation are well known so we merely state the final result here, which is quite close to the equations used in kivstamm in $1+1$ dimension, or the equations found in wen ; tsitsas ; walasik . Full derivation can be found in Appendix A. The equations of motion turn out to be:

[TABLE]

Here, $\nabla_{\perp}\equiv(\partial_{x},\partial_{y})$ is the nabla operator in the transverse plane, $k$ is the wavevector along the $z$-direction, $b$ is the characteristic propagation length along the $z$-axis. Equations of motion (4-5) together with the equations (1-3) for the permittivities contain the following five parameters: $\epsilon_{D0},F,\Gamma,\gamma,\omega_{0}$. Realistic values for all the parameters are discussed in sadriv05 . The natural length scale of the model is dominated by the $1/\omega_{0}$ scale. Dimensional analysis of the terms on the right-hand side of (4) determines the length scale $b$ in (4-5) as $b\sim 1/\omega_{0}$. Both in analytical and in numerical calculations, we express the transverse coordinates ($x,y$) in millimeters but the longitudinal coordinate $z$ is often stated in units of $b$. This is because the length scales of all patterns in the transverse plane are similar whereas the propagation lengths along $z$ can vary by an order of magnitude as $\gamma$ and $\Omega$ are varied, so it is more natural to express them in terms of the characteristic distance $b$.

II.1 A field-theoretical model

For some theoretical considerations it is useful to formulate a Lagrangian (gradient) model which captures the essential features of the equations of motion (4-5). As it often happens in studies of complex nonlinear pattern-forming systems, we cannot easily write the original equations in such a form. Instead, we construct a field theory which yields equations of motion somewhat different from the original ones but which still give the same phenomenology, and are able to explain the results of numerical calculations with the equations (4-5).

Let us think what such a field theory would look like. The terms with the gradient of magnetic permittivity obviously introduce dissipation, which physically originates from the losses in the inductive rings of the metamaterial. In general, dissipative systems do not have a Lagrangian, although a number of generalized Lagrangian approaches exist for dissipative systems: either with more general functional forms of the Lagrangian, or with a dissipative function in addition to the Lagrangian, or with extra degrees of freedom disslag1 ; disslag2 . We will take the first, most conventional of the three approaches: we will consider a conventional Lagrangian (no dissipative function, no extra degrees of freedom) which gives slightly generalized equations of motion compared to (4-5), with dissipative terms for both electric and magnetic fields coming from the complex terms in the effective potential. The effective action reads:

[TABLE]

The last term in $\mathcal{L}_{E}$ equals $-\omega^{2}\int_{0}^{EE^{\dagger}}dx\epsilon(x)$, analogously to the corresponding term in $\mathcal{L}_{H}$, but since $\epsilon$ is polynomial in $E^{\dagger}E$, the integral can be solved explicitly. Now (6) gives the equations of motion:

[TABLE]

where $\Phi_{E,H}$ are related to the fluxes of the electric and magnetic field (prime denotes the derivative of $\epsilon$ and $\mu$ with respect to their arguments $E^{\dagger}E$ and $H^{\dagger}H$):

[TABLE]

These are the extra terms compared to the physical equations (4-5).111The dissipative term proportional to $\nabla_{\perp}H$ in (8) is also absent in the original equations, but that one is easy to interpret: we make both $\mathcal{L}_{E}$ and $\mathcal{L}_{H}$ complex, so both fields have dissipative dynamics. Inserting $\partial_{z}E^{\pm}$ from the equations of motion (7-8) into the above we derive:

[TABLE]

and analogously for $\Phi_{H}$, with $\epsilon\leftrightarrow\mu,E\leftrightarrow H$. This term is proportional to a total derivative, and is therefore related to the flux $(E\nabla_{\perp}E^{\dagger}-E^{\dagger}\nabla_{\perp}E)/\mu$. For slowly-changing $\epsilon$ and $\mu$, which is often the case in our system (i.e., for $\epsilon^{\prime},\mu^{\prime}\ll\epsilon,\mu$), this term is small, which partly justifies the choice (6) for the Lagrangian. But the ultimate justification, as it frequently happens, is that a posteriori we will find that this model is able to explain the features observed in the numerics. Therefore we will not try to interpret the term (11) in detail.

III Geometry and stability of vortices

We will now sum up our numerical results which demonstrate the breaking of the circular symmetry of the vortex beams and their decay during the propagation. We always start from a Gaussian input beam with a topological charge $Q$, of the form $E(r,\phi;z=0)=E_{0}\times e^{-r^{2}/2\sigma^{2}}e^{iQ\phi}$ and analogously for the magnetic field, with amplitude $H_{0}$ but with the same vortex charge $Q$. Therefore, we always give an exact vortex as an input. The parameters of the model were chosen so that the permittivities $\epsilon$ and $\mu$, given by (1) and (2) respectively, are of order unity. This serves to limit the dissipative effects, so that the propagation along the longitudinal direction can be clearly observed. Same phenomena are found for arbitrary values of $\epsilon$ and $\mu$ but on a different length scale. We do not aim at a stamp-collecting exhaustive description of patterns for all possible parameter values, so we will focus on just a few relevant cases. We are mainly interested in left-handed materials ($\epsilon,\mu<0$) and how they compare to right-handed ones, so for the dielectric constant we always choose the self-defocusing Kerr nonlinearity ($\alpha=-1$) with a linear part $\epsilon_{D0}=12.8$, which has both a left-handed and a right-handed regime. To check the effects of dissipation, we either adopt $\gamma=0$ in (1), i.e., the lossless case, or we tune $\gamma$ so that $\omega_{0}^{2}/(\omega^{2}+\imath\gamma\omega)=1/2$. In other words, we impose either $\Im\epsilon=0$ or $\Im\epsilon=\Re\epsilon$. This is for simplicity and to avoid probing a huge parameter space for all possible $\gamma$ values; from now on we will call these cases simply lossless $\epsilon$ and dissipative $\epsilon$. The filling factor is $F=0.4$ and the magnetic dampening coefficient is $\Gamma=10^{9}$ Hz; these values are kept fixed in all calculations. Numerical calculations are performed with an operator split algorithm described in detail in the Appendices of nasstari .

The nonlinear frequency of the oscillator rings is obtained as a solution to (3). Of the three branches of the solution, we take the one that yields a negative real value of $\mu$ for $\omega>\omega_{0}$ (Fig. 1). We have freely taken $\omega=9.8\times 10^{9}$ Hz to represent a left-handed medium, and $\omega=7.0\times 10^{9}$ Hz to represent a right-handed medium. The transverse profiles are displayed in Fig. 1. We see there is a well-defined left-handed regime.

Now we discuss the transverse intensity profile for different initial beam configurations, with vortex input beams as explained in the beginning. We observe the following features:

Circular symmetry of the vortex input always breaks down to a discrete group.

1. (a)

In a dissipative left-handed medium, the discrete symmetry group for a vortex of charge $Q$ is $C_{3Q}$, before breaking down to simple $C_{2}$ axial symmetry at longer distances – Fig 2(a). 2. (b)

In a dissipative right-handed medium, the discrete symmetry group for a vortex of charge $Q$ is $C_{2Q}$, before breaking to $C_{Q}$ and then to $C_{2}$ axial symmetry at longer distances – Fig. 2(b). 3. (c)

In a lossless left-handed medium, the discrete symmetry group for a vortex of charge $Q$ is $C_{3Q}$ for very short distances, before quickly breaking down to $C_{Q}$ and finally $C_{2}$ – Fig. 2(c). 4. (d)

In a lossless right-handed medium, the discrete symmetry group for a vortex of charge $Q$ is $C_{2Q}$, before breaking to simple $C_{2}$ axial symmetry at longer distances – Fig. 2(d). 2. 2.

Vortices decay approximately exponentially as they propagate along the longitudinal axis. Fig. 4 shows the intensity of the beam across the $z$-axis, for various regimes. At early $z$ values, total intensity may behave non-monotonically and non-universally but on longer scales it decays exponentially. For different charges, the intensity curves collapse to a unique exponential function at large $z$. As could be expected, lossless and dissipative cases differ somewhat and collapse to different curves.

The bottom line is that there is a vocabulary of patterns with $C_{Q},C_{2Q},C_{3Q}$ symmetries. One of them dominates in each case (left/right handed, dissipative/lossless) but at longer propagation distances the symmetry can change, before the intensity drops to near-zero from dissipation. The final stadium of $C_{2}$ symmetry is only seen at very low intensities, so it might be practically unobservable in experiment; that is why we say the vocabulary only has three possible patterns, excluding $C_{2}$.

The findings above are further corroborated by Fig. 3 which shows the vortices with different charges $Q=1,2,3$ in the same regime (dissipative left-handed, (a), and lossless left-handed, (b)). As claimed above, the symmetry is $C_{3Q}$ in the panel (a), and (except at small $z$ values) $C_{Q}$ in the panel (b). Finally, it is obvious that there is some mixing of patterns: the polygons are never exactly regular, so the groups $C_{n}$ are certainly not exact symmetries; we use the $C_{n}$-nomenclature merely for convenience.

One interesting phenomenon in Fig. 2(c) is that the pattern rotates along the $z$ axis. This can be understood as excitation of multiple angular modes (of the form $e^{\imath lz}$ with various $l$-numbers) as the beam travels along the sample. This is a well-known consequence of nonlinear terms rmp ; book and typically depends on the relative strength of nonlinear mode interactions compared to energy density $|E|^{2}+|H|^{2}$ and dissipation $\gamma$. We will not explore it in quantitative detail in this paper as it is only tangential to our main topic of radial symmetry breaking; as one can see, the structure remains the same; just the orientation changes.

One might rightly worry that the initial conditions which contain a vortex in both electric and magnetic field are not very realistic, as in most materials the electric field dominates the optical response.222We thank an anonimous referee for pointing out this issue. Therefore, in experimental practice, one typically prepares a vortex in the electric field making use of phase masks or some other method, and the initial magnetic field distribution is completely analytic. In Appendix B we repeat the calculations from Fig. 3 and show that the outcome is the same, including the vocabulary of patterns and their $C_{n}$ shapes. Therefore, the $E$-$H$ symmetric ansatz is merely a matter of convenience, and the realistic regime where $|H|\ll|E|$ is in fact covered by our work.

Fig. 4 shows that at long times the decay of intensity is universal for given dielectric dampening coefficient $\gamma$, which suggests the main mechanism of dissipation is in fact the radiative loss. This is because we deliberately chose $\epsilon,\mu$ with small imaginary parts (for $\epsilon$ it can also be zero), so the losses in the medium are not so important when it comes to total energy (they are still important for being nonlinear and influencing the patterns). One important difference between the lossless medium (black and blue symbols in Fig. 4) and the dissipative medium (red, magenta) is that the former has a short interval of growing intensity, before reaching the universal regime of radiative decay. The physical reason is that the polarization, i.e., the rearrangement of charges in the self-defocusing metamaterial reduces the overall electrostatic potential energy of the medium, and this energy becomes available to the beam, increasing its intensity. Clearly, once the radiative losses overcome the total potential energy available, the intensity decays. The growth is clearly a transient effect which cannot persist for long $z$-intervals. A formal way to understand this is that the nonconservation of energy is encoded by the last term in (4), which can have a positive or negative imaginary part depending on the sign of $\partial_{z}\mu/\mu$. At large values of $z$, we expect to enter a universal regime where this sign is constant, because the radiation loss dominates over nonlinearities and the exchange of energy between the beam and the medium; this is the universal decay regime in the figure.

IV The theory of vortex evolution

The phenomenology described in the previous section can be understood on several levels. At the crudest level, we can introduce a variable-separation ansatz in the equations of motion and then linearize them in the amplitude (but not in the phase!). This picture explains the $C_{2Q}$ patterns, but not the $C_{3Q}$ and $C_{Q}$ regimes. It also does not explain the instabilities, that is the changes and disappearance of patterns during the $z$-propagation. For the full picture it is necessary to take into account the nonlinear effects through the loop corrections, i.e., to move perturbatively beyond the amplitude-linearized solution. A qualitative insight of the symmetry breaking can however be obtained also in a simpler and more elegant way, directly from the symmetry analysis of the model Lagrangian (6). Therefore, after finishing the amplitude-linearized analysis and the loop corrections from nonlinearity, we will obtain the same results from a unified mean-field treatment of the (nonlinear) model Lagrangian.

A note on terminology is in order. The solutions we find are not the textbook type of vortex with phase dependence solely of the type $e^{\imath Q\phi}$; rather, the dependence on the phase is more complicated, i.e., the phase is doing more than just the winding, but it is still true that the circulation of the phase around some point (the location of the vortex core) is an integer – the topological charge of the vortex. Such solutions are sometimes called spirals rabin whereas the term vortex is reserved for the simple winding-phase solutions. We nevertheless stick to the widespread term ”vortex” for any topologically charged solution under the fundamental group of the $U(1)$ phase symmetry.

IV.1 Amplitude-linearized solution

We will separate variables in the equations of motion (4-5) (or the Lagrangian equations (7-8), which do not differ from the original equations at the amplitude-linearized level) and then plug in the vortex ansatz. The vortex ansatz is a solution which has a winding phase $\Phi$ with some winding number $Q$, for a constant (averaged) value of the permittivity $\mu_{c}=\mathrm{const.}$, because we ignore the nonlinear dependence of $\mu$ on $|H|$. The vortex solution of winding number (topological charge) $Q$ in cylindrical coordinates $(r,\phi,z)$ can be separated into regular and vortex parts:

[TABLE]

We represent the vortex part as

[TABLE]

and analogously for the magnetic field. Along the $z$-axis we get $Z_{E}(z)=e^{\imath\lambda z}$ as expected, and the eigenvalue $\lambda$ is arbitrary for now, i.e., it is determined by the boundary conditions along the $z$-axis. Upon inserting (13) into (4), the equation separates into the angular part and the radial part. The former reads

[TABLE]

where $l$ is the eigenvalue of the angular part. This is the crucial equation – the phase dynamics is nonlinear because $\mu$ is in general complex and the terms with $\nabla_{\perp}\mu$ contain nonlinear dependence on the phase. The equation is easily solved by first introducing $w\equiv\Phi^{\prime}$ and then reducing it to quadratures. The outcome is

[TABLE]

In other words, we still stay with a winding solution but various winding numbers (equal to $\sqrt{Q^{2}+l^{2}}$) are possible when multiple modes are excited. Clearly, only the solutions with integer windings are physical, otherwise they would not be single-valued. The most general solution is thus a superposition of solutions $Z_{E}(\lambda,l;z)\Phi_{E}(\lambda,l;\phi)R_{E}(\lambda,l;r)$ with different $l$-modes so as to result in a single-valued function. Now the radial part acquires the form

[TABLE]

with $\tilde{\omega}\equiv\omega(1-\omega_{0}^{2}/(\omega^{2}+\imath\omega\gamma))$. If we disregard the cubic term (amplitude-linearized approximation!),333This is justified at least in some interval of $z$ values, as the system is dissipative and loses power $\int(E^{2}+H^{2})$, so the amplitude progressively decreases along $z$. the well-known solution in terms of Bessel functions is obtained:

[TABLE]

Here, $J$ and $Y$ are the Bessel functions of first and second kind, respectively. Similar solutions $Z_{H}(z),\Phi(\phi),R_{H}(r)$ are obtained for the magnetic field. The angular equation is identical for both fields: for this reason we have one solution $\Phi$ for both $E$ and $H$. The eigenvalues $\lambda,l$ and the values of the constants $c^{(1,2)}_{E,H}$ are determined by the boundary conditions. Obviously, (15) imposes the $C_{2Q}$ symmetry, if $l=0$. This simplest case is not necessarily the stable solution. We might have a sum over many $l$-values. In principle, such sums may yield more complicated patterns, however we will see that when the physically reasonable boundary conditions are implemented (decay at infinity, single-valuedness everywhere) one typically always has the robust $C_{2Q}$ pattern. One important consequence of the fact that multiple $l$-modes are possible is that due to nonlinear effects a new $l$-mode can be created during the propagation along the $z$-axis. We have already seen an example in Fig. 2(c). A quantitative analysis of this phenomenon requires a full nonlinear model and so can only be studied within the formalism of the next section.

This solution is not very satisfying but reproduces some of the features from the numerics, summarized at the start of the previous section: (1) the reduction of the full $O(2)$ symmetry down to a discrete symmetry $C_{n}$ for some $n\in\mathbb{N}$, i.e., the polygonal form of the vortex (2) the value $n=2Q$ is true in some but not in all situations. We show the solutions for a single angular mode from (15,17) in Fig. 5(a). In Fig. 5(b), we show a linear combination of angular modes with $l=0,1,2$, with the coefficients $c^{(1,2)}_{E,H}$ in (17) chosen so that the total intensity still decays sufficiently fast at infinity. The symmetry is still $C_{2Q}$. Apparently, the regimes with the $C_{Q}$ and $C_{3Q}$ symmetries require loop corrections from nonlinear $\mu$ to be taken into account.

IV.2 Loop corrections

The origin of the breaking of radial symmetry is the fact that a discrete set of modes in Fourier space is selected. This is best seen from the Fourier transform of the solutions (15) and (17). We will calculate the propagator $G(\mathbf{u})$ at constant $z$, i.e., the Fourier transform $\mathbf{r}\mapsto\mathbf{u}$ of the solution with a Dirac delta source. This source imposes the boundary condition $R_{E}(0)\to\infty,\int drr\cos\phi R_{E}(r)=1$, giving $c^{(1)}_{E}=0,c^{(2)}_{E}=2\pi/\Gamma(Q/2)$ in (17). Fourier-transforming $(x,y)\mapsto(u_{x},u_{y})$ we get for a single mode (17), making use of the Bessel and Lommel integrals:

[TABLE]

Here, $\Lambda$ is the ultraviolet (UV) small-length/high-momentum cutoff, i.e., the Fourier transform is performed by integrating $\int_{1/\Lambda}^{\infty}dr\int_{0}^{2\pi}d\phi$. The cutoff has a clear physical meaning: $1/\Lambda$ is the size of the vortex core (where the vortex ansatz stops working because the gradient of the field becomes too high). We clearly do not get anything new by just Fourier-transforming. The goal is to move beyond the amplitude-linearized approximation of the previous section by considering the effects of non-constant permittivity $\mu$ instead of constant (averaged) $\mu_{c}$. This calculation is essentially elementary but might be tedious and boring for readers who are not fond of perturbative field theory. Most of the integrations are in Appendix C. Even the rest of this subsection can be skipped until the the equation (22) where we discuss the final result.

Putting $\mu$ from (3) in place of $\mu_{c}$ requires the solutions for $\omega_{0NL}$ in terms of the magnetic field. The solutions are readily found from the Cardan formulas (we do not give them explicitly as they are cumbersome and not very illustrative). But the form of the $H$-dependence of $\omega_{0NL}$ is seen already from the Viete formula:

[TABLE]

so the solutions depend on $|H|^{2}$ only, with no higher powers of the magnetic field. Inserting this into $\mathcal{L}$, we get the nonlinear correction of the form:

[TABLE]

We thus have two quartic interaction terms and two quadratic terms. We do not intend to calculate the loop corrections in full detail; it is not worth the effort as we only want to capture the symmetry, i.e., the form of the angular dependence. First of all, the quadratic corrections $g_{2,0,0},g_{0,2,0}$ trivially renormalize the parameters in the bare propagator and do not change its functional form. Lowest-order nontrivial loop corrections to the self-energy come from $g_{0,2,2}$ and $g_{2,0,2}$. The electric field receives the correction $G_{E}^{-1}\mapsto\left(G_{E}+\Sigma_{E}^{(1)}+\Sigma_{E}^{(2)}\right)^{-1}$ with

[TABLE]

We will write all equations for $E$, because this field receives interesting corrections from the gradient of $\mu$ in Eqs. (4,7). The magnetic field does not couple to the permeability $\epsilon$ in the same way in the original equation (5), and in the Lagrangian form (8) it does but $\epsilon$ does not contain such strong (non-polynomial) nonlinearities as $\mu$. One- and two-loop corrections appear not only in the self-energy but also in the vertex operators. However, the vertex corrections only have a weak momentum dependence and consequently the coordinate dependence (geometric patterns) of the solution is not significantly affected by them. For that reason we will not discuss them in detail.

The correction $\Sigma_{E}^{(1)}$ is the Hartree correction with a single vacuum bubble which is not very interesting: it merely introduces an additional mass term and does not influence the momentum dependence and thus the geometry of the patterns. As could be expected from power counting, it is logarithmically divergent in the UV cutoff $\Lambda$. Of course, this is not a problem in an effective theory; we have already explained the physical meaning of $\Lambda$. The watermelon diagram $\Sigma_{E,H}^{(2)}$ is crucial: it is momentum-dependent. Its calculation is found in Appendix C. An informal way to estimate its effect is the following: the leading contribution comes from the region where $\mathbf{u}\approx\mathbf{u}^{\prime}-\mathbf{u}^{\prime\prime}$ because this is a pole of the self-energy correction. Then we are left with angular integrals only, and they reduce to integrals of products of three rational functions (for the three propagators in (21)) of the half-angle – this gives rise to $3\phi/2$ in the argument of the cosine. Now the dressed propagator $\left(G_{E,H}^{-1}+\Sigma\right)^{-1}$ needs to be Fourier-transformed back to real space. We will only do this approximately (it is likely impossible to do exactly in closed form). The outcome is

[TABLE]

No doubt the reader sees that the terms $\cos(3Q\phi/2),\sin(3Q\phi/2)$ give a pattern $|E_{\mathrm{vort}}|^{2}$ with $3Q$ branches, in addition to the $2Q$-polygons obtained from the term $\cos(Q\phi)$. The interference between the two patterns might (1) break the symmetry completely (2) lead to $C_{Q}$ symmetry if the relative phase between the leading term and the corrections is approximately $2\pi/Q$. Both cases are seen in numerical work: $C_{3Q}$ appears in all left-handed materials (Fig. 2(a,c)), and elements of $C_{Q}$ symmetry are present in almost all cases at long propagation distances $z$ (Fig. 2(a,b,c), Fig. 3).

The self-energy has an imaginary part (equivalently, the solution (22) exhibits exponential decay in $z$), meaning that these configurations are not stable - they are only seen up to some propagation distance $z$. The exact order (along $z$) and stability of each of the patterns depends on the details of the permeability $\epsilon$. One important and universal lesson is however that the decay rate (the real part of the exponent in (22)) is proportional to $Q^{2}$, therefore the higher the value of $|Q|$, the faster it decays. This supports the general intuition that vortices with high winding numbers are not stable. But unlike the simplest case of the XY model or a superfluid where the stability only allows $Q=\pm 1$, we can in principle have arbitrarily high $Q$ as we have seen also in the numerics; their lifetimes are smaller and smaller as $Q$ grows, but still finite. The exponential decay itself is also confirmed by the numerics, as seen from Fig. 4.

IV.3 Isotropy breaking – the look from the action

The basic mechanism leading to the symmetry breaking $O(2)\mapsto C_{3Q}\mapsto C_{2Q}\mapsto C_{Q}$ is seen already from the model Lagrangian (6). The symmetry breaking is essentially the consequence of the interplay of the nonlinear-sigma-model form of the kinetic term and the complex nonlinearity of the magnetic permittivity $\mu$. Therefore, we can take a static approximation of the $z$-dynamics, ignoring the $z$-dependence; clearly, in that framework we can only obtain the vocabulary of patterns, not the relative stability of $C_{Q},C_{2Q},C_{3Q}$.444We could take the ansatz $e^{\imath\lambda z}$ instead; it would merely modify $k^{2}\mapsto k^{2}-\lambda$. The separation of variables remains a natural ansatz, and the vortex nature of the solution implies $E_{\mathrm{vort}}=E_{0}(r)e^{\imath\Theta(\phi)}$ with $\oint d\phi\Theta(\phi)=2\pi Q$ and analogously for the magnetic field. The Lagrangian (6) then becomes:

[TABLE]

The fact that $\mu$ contains $\omega_{0NL}^{2}(|H|^{2})$, which is in turn the solution of the cubic equation, introduces a branch cut in $H$ because of the cubic roots. This is the simplest explanation of the origin of the $C_{3Q}$ symmetry. More quantitatively, the story follows exactly the Landau-Ginzburg paradigm: while the initial Lagrangian only depends on $|E|^{2}$ and $|H|^{2}$ and thus preserves isotropy, the saddle-point solution is given by the equation

[TABLE]

where we have used that $\mu=\mu(\omega_{0NL}^{2})$ and $\omega_{0NL}^{2}=\omega_{0NL}^{2}(|H|^{2/3},|H|^{4/3})$ (from the Cardan formulas). With the ansatz adopted above, the amplitude equation for $E_{0}(r)$ is the nonlinear amplitude equation (16). The equation for the phase part $\Theta$ is more interesting. It reads

[TABLE]

The cubic root carries a branch cut, and the last term really evaluates to $2\mu^{\prime}/3\mu^{2}\times H_{0}^{-1/3}e^{-\imath\Theta/3+2n\pi\imath/3}$ with $n=-1,0,1$. The solution $\Theta_{0}$ which satisfies the phase winding condition is obtained in implicit form as

[TABLE]

where $K_{n}$ is a constant determined by the amplitude solution and depending also on $n=-1,0,1$; its exact value is hard to find analytically as we do not know the solution to the amplitude equation in the nonlinear regime. But that is not crucial for our general argument. The point is that the system can choose a solution with any of the values $n=-1,0,1$, i.e., even though the equations of motion (and the Lagrangian) are isotropic, the solution is not. Each $n$-branch behaves as $\sim 1/\cos^{2}(Q\phi/2)$, only they are rotated by $\pm 2\pi/3$ with respect to each other; and each of them has a $C_{Q}$ symmetry. Put together, the three branches give $C_{3Q}$ patterns. But all that holds if two of the cubic roots are complex. If all cubic roots are real, the phase remains single-valued, and we only have $C_{Q}$ symmetry, coming directly from (26) if we fix $n=0$, i.e., if we only keep a single branch.555We use the fact that a cubic equation has either one or all three solutions real.

What is the regime in which cubic roots are real and the symmetry is $C_{Q}$, as opposed to the complex roots and $C_{3Q}$ patterns? The easiest way is to look at the cubic equation (3) for the magnetic permeability (and the nonlinear frequency $\omega_{0NL}$). For $\mu>0$ (right-handed regime), the roots are all real and $C_{3Q}$ patterns cannot occur. Indeed, the $C_{3Q}$ phase is only present in Fig. 2(a,c), in left-handed media.

This approach is much more physical and elegant than the tour-de-force calculations of the previous two sections but it does not give explicit solutions for $E,H$; it only classifies the symmetries of the solution. This is why we we still needed the perturbative linear and two-loop analysis, to arrive at more quantitative results.

The saddle-point solution (26) is nonlinear, unlike the linearized solution found in the first subsection (15). It is not a vacuum in the usual field-theory sense however, as it is not constant. We are dealing with dynamical criticality of the kind discussed in cross . In the vicinity of this solution, the Lagrangian describes the fluctuations of amplitude $\delta E,\delta H$, and the fluctations of phase $\delta\Phi$. Similar to the $O(3)$-type spin models vortbook and multi-beam optical systems nasstari , and unlike simple XY-type models, the phase and amplitude fluctuations mix. By analyzing the fluctuation equations, it should be possible to understand analytically also the transition from the left-handed to the right-handed regime as the parameters are varied, i.e., what are the instabilities that drive it. We will not attempt that here; it is a long subject that deserves separate work.

V Toward experimental verification and applications

We will now briefly discuss what an experimentalist can learn from our results and what to look for in practical work. Wave propagation through the metamaterial can be observed by measuring the transmission coefficients $S_{ij}$. From these coefficients, one can also reconstruct the electric field intensity $|E|^{2}$, which can be directly compared to our intensity maps like Figs. 2, 3 alitalo . Another quantity which can be measured is the voltage waveform, which can be used to construct amplitude envelopes kozyrev .

Therefore, the predicted symmetry breaking is in principle directly observable. But the question remains how widespread it will be for realistic values of the parameters. From a more applied viewpoint, this question is reversed: how to make a vortex transmission through a left-handed waveguide stable. In other words, how not to observe the symmetry breaking. It is true that the phenomenon disappears as soon as the vortex charge is zero, i.e., when the beam is not a vortex. However, the vortex patterns are likely important in applications. First, as a topologically protected object with conserved charge, a vortex is among the natural candidates for computational devices and information transmission (for the same reasons that solitons are also interesting in that regard: they are robust to noise, carry a discrete ”quantum” number, i.e., charge and are stable to small local perturbations). Second, in the presence of impurities in the sample, vortices can form in a nonlinear metamaterial from the initially non-vortexing beam vortbook .

Let us focus on the left-handed regime, which is the most interesting and the most relevant for applications. The first condition is therefore to be in the frequency regime with $\mu(\omega)<0$. This can be checked directly from Eq. (2) as we did in Fig. 1(c). The second issue is that the symmetry breaking takes some finite time, i.e., some finite propagation length, which is of order $b$; as can be seen from Fig. 3 and directly from Eqs. (4-5), this is the length scale over which the patterns change. On the other hand, the one-loop calculation (22) shows that the intensity decays with the rate $\sim a^{-2Q^{2}}$. As long as this is less than the characteristic length $b$, one will likely not see the symmetry breaking but just eventual dissipation of the beam. Therefore, these two scales should be compared for some reasonable parameter values. We show this in Fig. 6(a) for $F=0.4$, $\epsilon_{D0}=12.8$, $\gamma=1$ GHz, $\omega_{0}=10$ GHz. Apparently, the length scale of the $C_{n}$ pattern development (red dotted line) is nearly always shorter than the dissipation scale (blue dashed line), so we expect that the effect predicted in the paper is readily seen in experiment, at least for $Q=\pm 1$. For larger vortex charges, the dissipation grows quickly and high $Q$ values are probably not easily observed. Conversely, if the goal is to keep a stable radially symmetric vortex pattern, one should remain at small frequencies, although for $\omega\ll\omega_{0}$ the material is not strongly left-handed, as can be seen from the $-\mu(\omega)$ dependence, also given in the figure.

There is still one remaining issue. Our theoretical approach, based on a pair of nonlinear Schrödinger-like equations, inherently disregards some effects. It describes a quasimonochromatic wave without wave mixing or dissipation due to higher harmonic generation rmp . Such phenomena become significant for strong nonlinearities, so we should compare the nonlinearities in $\epsilon$ and $\mu$ to the typical energy (frequency) scale of the vortex. In Eqs. (1-2) the approximate ratios of the nonlinear to linear terms are given by $|E|^{2}/\epsilon_{D0}$ and $\omega_{0NL}/\omega_{0}\sim(A/H)^{1/3}$. The first scale is frequency-independent and solely depends on the beam intensity. The second scale depends on frequency and needs to be inspected more closely. In Fig. 6(b) we plot the nonlinearity ratio for the magnetic field for a range of frequencies $\omega$, again together with the permittivity to make sure we are at the same time in the left-handed regime. The relative nonlinearity strength quickly saturates around a value $0.06\ll 1$, so we are rather confident that our equations of motion still make sense.

Altogether, the conclusion is that the breaking of radial symmetry is observable by standard means (measuring the transport coefficients and reconstructing the intensity map at the exit face of the metamaterial), as long as the frequency of the wave is not too low. This kind of instability kicks in at shorter propagation lengths (of order $0.1$ mm in Fig. 6(a)) than the nonlinear diffraction effects studied for breathers in boardman , suggesting that vortex signals are more fragile and less convenient for information transmission.

VI Discussion and conclusions

Our main result is contained already in the title – left-handedness and nonlinearity together create the breaking of the $O(2)$ symmetry down to a discrete group, with the pattern vocabulary consisting of the $C_{3Q},C_{2Q},C_{Q}$ patterns. In the right-handed system with the same nonlinearity the isotropy is broken again, but the pattern vocabulary only has $C_{2Q}$ and $C_{Q}$ stages. How exactly the patterns evolve into each other and through which instabilities is not universal, and it depends on the exact form of $\epsilon$ and $\mu$. In our model, the $\epsilon$-dependence is mainly encapsulated in the dissipation $\gamma$: the left-handed non-dissipative case is usually dominated by $C_{Q}$ after a much shorter $C_{3Q}$ phase, whereas the dissipative left-handed metamaterials most prominently show $C_{3Q}$ patterns. For the right-handed materials, non-dissipative and dissipative dynamics shows mainly $C_{2Q}$ and $C_{Q}$ patterns, respectively.

A detailed account of the pattern dynamics was only possible through numerical work. But the vocabulary itself – the existence of symmetries $C_{3Q},C_{2Q},C_{Q}$ – we were able to understand analytically. The dynamic Landau-Ginzburg picture reveals this as a consequence of the cubic root nonlinearity in the magnetic permittivity, and the fact that the cubic equation has either two complex roots in the left-handed regime, or all three real roots in the right-handed regime, and the presence or absence of dissipation in the electric permeability. In the framework of our field theory model, the second derivative of the free energy (on-shell Lagrangian, Landau-Ginzburg functional) likely has a jump when the symmetry changes. This is a strong encouragement that the phenomena we observe here, and in general the walk through the pattern vocabulary, can be understood from the viewpoint of order/disorder transitions.

Similar phenomena were studied also in fan ; kevrek and above all walasik , where $C_{3Q}$ necklaces were found, within a model of left-handed metamaterials given in fan and similar to ours. Clearly, we have not exhausted this subject; more research is still needed to fully understand the transition between different symmetries and their instabilities. Vortices in metamaterials seem to be a promising arena, as in a metamaterial the nonlinearity and the frequency band where the materials is left-handed can to some extent be tuned at will. Therefore, the phase diagram of collective vortex interactions can also be studied, and is an obvious topic for future work.

Acknowledgments

This work has made use of the Sci-Hub service. We are grateful to Milan Petrović and Mariya Medvedyeva for helpful remarks. We also thank the referees for some important and stimulating questions. Work at the Institute of Physics is funded by Ministry of Education, Science and Technological Development, under grant no. OI171017.

Appendix A Derivation of the equations of motion from the Maxwell equations

Start from the definitions $\hat{D}=\epsilon\hat{E},\hat{B}=\mu\hat{H}$ and the Maxwell equations in absence of external charges and currents ($\rho=\hat{j}=0$):

[TABLE]

We make the following assumptions:

Small gradients of the permittivities $\epsilon$ and $\mu$, so their second and higher derivatives are disregarded. Since $\omega\propto k$, it means that mixed derivatives of the from $\partial_{t}\nabla\epsilon$ are also disregarded. In other words, the characteristic length scale $l$ along the $z$-axis on which $\epsilon,\mu$ change is assumed to be large compared to the characteristic scale $b$ of the changes in $E,H$. 2. 2.

The time dependence is harmonic so $\partial_{t}=-\imath\omega$.

Acting on the last equation by $\nabla\times$ and making use of the identity $\nabla\times\nabla\times\hat{H}=-\nabla^{2}\hat{H}+\nabla(\nabla\cdot\hat{H})$, one gets for the left-hand side:

[TABLE]

where we used $\nabla\cdot\hat{B}=0$ and disregarded the second derivative of $\mu$. The right-hand side yields

[TABLE]

so we obtain

[TABLE]

For the $\hat{E}$ field we start from the third Maxwell equation, act by $\nabla\times$ and find for the left-hand side:

[TABLE]

and for the right-hand side we get

[TABLE]

so

[TABLE]

For our geometry we take the paraxial beam approximation, with the ansatz $\hat{E}=E(x,y)e^{\imath(kz-\omega t)},\hat{H}=H(x,y)e^{\imath(kz-\omega t)}$, so the nabla acts as

[TABLE]

and the Laplacian operator acts as

[TABLE]

and analogously for the magnetic field. Now to write the equations motion in the final from we rescale $E\to E\cdot 2kb,H\to H\cdot 2kb,z\mapsto z\cdot 2kb$, where $b$ is some characteristic length scale along the $z$-axis, and divide the equations by $bk^{2}$ to obtain the equations (4-5), reprinted here for convenience:

[TABLE]

For comparison to the equations given in zarov03 ; sadriv03 ; sadriv05 , one needs (1) to rescale $H\mapsto\omega^{2}/c^{2}H$ to get the term $-\gamma^{2}H=-k^{2}/\omega^{2}$ in (37) and (2) to absorb the factor $-k^{2}$ in (36) in the definition of $\epsilon_{D0}$. This is possible as $\epsilon$ and $\mu$ have a constant term (equal $\epsilon_{D0}$ and $1$ respectively) so the product $\epsilon\mu$ also has a constant term proportional to $\epsilon_{D0}$, and the contribution $k^{2}E$ can be absorbed as $\epsilon_{D0}\mapsto\epsilon_{D0}-k^{2}$. We thus arrive at a system identical to that from zarov03 , except for the extra terms for the propagation along the $z$-axis.

Appendix B Configurations with no vorticity in the magnetic field

Here we show that our results stay valid also when only the electric field has vortex patterns whereas the magnetic field starts analytic everywhere. As we discuss in the main text, this situation is experimentally more relevant than the one assumed in most calculations in the paper (that both the electric and the magnetic field have a vortex as they enter the material). The electric field is typically a few orders of magnitude more intense than the magnetic field, as seen in zarov03 . Therefore, one typically controls the electric field directly, imposing a given boundary condition at the front end of the material. Despite this fact, the magnetic field remains very important: the coupled equations of motion (4-5) require both $E$ and $H$ to be nonzero. Indeed, as explained in zarov03 , the left-handedness comes as a consequence of the hysteresis-type dependence of the magnetic permittivity on $H$. So while it is crucial that $E$ and $H$ are both nonzero, it is also true that the results should remain valid for $|H|\ll|E|$, and for the boundary condition that only has a vortex in $E$ at the front of the metamaterial, not for $H$. With such boundary conditions and the same parameter values as before, Fig. 7 repeats the calculations of Fig. 3. Obviously, the symmetries remain the same and the similarity of the results for the two cases is striking. Obviously, the $|E|^{2}$ map is insensitive to the details of the initial magnetic field pattern, as one expects from experiments and common wisdom in nonlinear optics. We are thus content that the numerically simplifying assumption of identical $z=0$ boundary conditions for $E$ and $H$ does not put into question the findings of our work.

Appendix C The calculation of the self-energy diagrams

We discuss here in some more detail the equations (21) from the main text. First we give the expressions for the couplings $g_{2,0,0},g_{0,2,0},g_{2,0,2},g_{0,2,2}$, which come from the expansion over the magnetic field $H$ of the nonlinear dependence $\mu(H)$ in (20):

[TABLE]

For simplicity, we will treat the case when $\lambda=k$ and thus $g_{0,2,0}=g_{0,2,2}=0$. This simplifies the calculations substantially while it does not change the symmetry of the solution. It is possible to evaluate the diagram $\Sigma^{(1)}$ exactly in terms of sine and cosine integrals $\mathrm{Si}$ and $\mathrm{Ci}$. The angular integration is straightforward, the integration over $u$ results in four combinations of the trigonometric integrals, for the four terms in (18). Three of the four integrals are finite and therefore they just shift the mass term. The third term of the propagator is logarithmically divergent:

[TABLE]

To judge the effect of this term, we should extract the mass squared $r_{m}$ of the bare propagator, writing it out for small $u$:

[TABLE]

Since $G^{-1}(\mathbf{u}\to 0)\propto u=0$, the bare propagator is massless. The one-loop correction $\Sigma^{(1)}$ therefore gives a cutoff-dependent mass $r_{M}\sim\log\Lambda$, which could be absorbed in the overall normalization of the propagator. As we declared in the main text, the one-loop self-energy does not do much.

The crucial diagram $\Sigma^{(2)}$, the popular watermelon diagram, cannot be calculated exactly. It can be evaluated in the regime of small external momentum $\mathbf{u}$, i.e, when $u<u^{\prime},u^{\prime\prime}$; more precisely, we can look at the regime when $u<u_{0}<u^{\prime},u^{\prime\prime}$ for some (arbitrary) scale $u_{0}$ and expand in a series in $u/u_{0}$. Let us denote such entity by $\Sigma^{(2)}(\mathbf{u};u_{0})$: it contains enough information for our purposes: we are interested mainly in angular integrations which determine the symmetry, and these can be done exactly as they separate from the integrations over the module $u$ in the small-$u$ limit. For $\mathbf{u}=0$ the watermelon diagram reads (with $\int\equiv\int_{0}^{2\pi}d\phi^{\prime}\int_{0}^{2\pi}d\phi^{\prime\prime}\int du^{\prime}\int du^{\prime\prime}$):

[TABLE]

One angular integration is performed by taking $\phi^{\prime}\mapsto\phi^{\prime}+\phi^{\prime\prime}$, which makes the $\phi^{\prime\prime}$ integral completely trivial, and the $\phi^{\prime}$ integral is evaluated in terms of the elliptic integrals $E,K$. The outcome is finite, hence it is observable (not only at the cutoff scale), and reads:

[TABLE]

In particular, this means that a nontrivial mass term is acquired, of the order $a^{3/2}$. This mass is anisotropic, and the factor $\cos(3Q\pi/2)^{2}$ is all we need for the $3Q$-polygon. The leading correction in $u/u_{0}$ is in fact inessential for the symmetry, but it is important as it contains a nonzero imaginary part, introducing a finite lifetime for such patterns. It reads

[TABLE]

At leading order, this tedious expression behaves like $1/r^{3}$, falling off much quicker than the bare propagator (18), which goes as $1/\sqrt{r}$ (most obvious from the Bessel-function form of the real-space solution), suggesting that the shape of the vortex, which is mainly determined by long-distance behavior, is not much influenced by the finite-$u$ correction to $\Sigma^{(2)}$.