Quantum criticality in photorefractive optics: vortices in laser beams and antiferromagnets

TL;DR

This paper explores vortex patterns in nonlinear photorefractive optics, revealing phase transitions analogous to condensed matter systems, including novel phases like spin-glass-like states, using statistical field theory methods.

Contribution

It introduces a field-theoretical analysis of vortex phases in photorefractive optics, connecting optical pattern formation to condensed matter models such as doped antiferromagnets.

Findings

Identification of multiple vortex phases including perfect conductor and frustrated insulator.

Discovery of a spin-glass-like phase in disordered photorefractive systems.

Mapping of optical vortex behavior to doped O(3) antiferromagnet model.

Abstract

We study vortex patterns in a prototype nonlinear optical system: counterpropagating laser beams in a photorefractive crystal, with or without the background photonic lattice. The vortices are effectively planar and described by the winding number and the "flavor" index, stemming from the fact that we have two parallel beams propagating in opposite directions. The problem is amenable to the methods of statistical field theory and generalizes the Berezinsky-Kosterlitz-Thouless transition of the XY model to the "two-flavor" case. In addition to the familiar conductor and insulator phases, we also have the perfect conductor (vortex proliferation in both beams/"flavors") and the frustrated insulator (energy costs of vortex proliferation and vortex annihilation balance each other). In the presence of disorder in the background lattice, a novel phase appears which shows long-range correlations and absence of long-range order, thus being analogous to spin glasses. An important benefit of this approach is that qualitative behavior of patterns can be known without intensive numerical work over large areas of the parameter space. More generally, we would like to draw attention to connections between the (classical) pattern-forming systems in photorefractive optics and the methods of (quantum) condensed matter and field theory: on one hand, we use the field-theoretical methods (renormalization group, replica formalism) to analyze the patterns; on the other hand, the observed phases are analogous to those seen in magnetic systems, and make photorefractive optics a fruitful testing ground for condensed matter systems. As an example, we map our system to a doped $O(3)$ antiferromagnet with $\mathbb{Z}_2$ defects, which has the same structure of the phase diagram.

Figures

Equations

$$\pm ı{\partial }_z{\Psi }_{\pm }(z;x,y;t)+\Delta {\Psi }_{\pm }(z;x,y;t)=\Gamma E(z;x,y;t){\Psi }_{\pm }(z;x,y;t),$$

$$\frac{\tau }{1+I(z;x,y;t)}{\partial }_{t}E(z;x,y;t)+E(z;x,y;t)=-\frac{I(z;x,y;t)}{1+I(z;x,y;t)}.$$

$$E(z;x,y;u)=-\frac{{\Psi }^{†}\Psi +I_x-\tau E_0}{1+\tau u+I_x+{\Psi }^{†}\Psi }=-1+\frac{1+\tau u+\tau E_0}{1+\tau u+I_x+{\Psi }^{†}\Psi }.$$

$$\mathcal{L}=ı{\Psi }^{†}{\sigma }_3{\partial }_z\Psi -|\nabla \Psi {|}^2+\Gamma {\Psi }^{†}\Psi -\Gamma (1+\tau E_0+\tau u)\log (1+\tau u+I_x+{\Psi }^{†}\Psi ),$$

$$V_{\mathrm{eff}}({\Psi }^{†},\Psi )=-\Gamma \log \frac{e^{{\Psi }^{†}\Psi }}{{\left(1+\tau u+I_x+{\Psi }^{†}\Psi \right)}^{1+\tau (E_0+u)}},$$

$${\Psi }_{0\pm }(\mathbf{r})={\psi }_{0\pm }(r)e^{ı\delta {\theta }_{\pm }(\varphi )+ı{\theta }_{0\pm }(\varphi )}.$$

$${\mathcal{H}}_{\mathrm{vort}}=\sum _{i<j}\left(g{\vec{Q}}_i\cdot {\vec{Q}}_j+g^{\prime }{\vec{Q}}_i\times {\vec{Q}}_j\right)\log r_{ij}+\sum _i\left(g_0{\vec{Q}}_i\cdot {\vec{Q}}_i+g_1{\vec{Q}}_i\times {\vec{Q}}_i\right).$$

$$\mathcal{Z}=\sum _{Q_{+},Q_{-}}{\left(\frac{\Lambda }{a}\right)}^2{e}^{-L{\mathcal{H}}_1}=\sum _{Q_{+},Q_{-}}{e}^{2\log \frac{\Lambda }{a}-L\left(g\vec{Q}\cdot \vec{Q}+g^{\prime }\vec{Q}\times \vec{Q}\right)\log \frac{\Lambda }{a}}$$

$${\mathcal{F}}_0^{(1)}=\left(g-\frac{2}{L}\right)\log \frac{\Lambda }{a},\vec{Q}=(\pm 1,0)\text{or}\vec{Q}=(0,\pm 1)$$

$${\mathcal{F}}_1^{(1)}=\left(2g-g^{\prime }-\frac{2}{L}\right)\log \frac{\Lambda }{a},(Q_{+},Q_{-})=(\pm 1,\mp 1)$$

$${\mathcal{F}}_2^{(1)}=\left(2g+g^{\prime }-\frac{2}{L}\right)\log \frac{\Lambda }{a},(Q_{+},Q_{-})=(\pm 1,\pm 1).$$

$$\frac{\delta \mathcal{Z}}{\mathcal{Z}}=1+\frac{y^4}{4}\sum _{q_{\pm }}\int d{r}_{12}{r}_{12}^3{e}^{g\vec{q}\cdot \vec{q}+g^{\prime }\vec{q}\times \vec{q}}{\left[\int dr{r}^2\left(g{\vec{Q}}_1\cdot \vec{q}+g^{\prime }{\vec{Q}}_1\times \vec{q}\right)\nabla \log |{\mathbf{R}}_1-\mathbf{r}|+\left(g{\vec{Q}}_2\cdot \vec{q}+g^{\prime }{\vec{Q}}_2\times \vec{q}\right)\nabla \log |{\mathbf{R}}_2-\mathbf{r}|\right]}^2.$$

$$\frac{\delta \mathcal{Z}}{\mathcal{Z}}=1+y^4\left(8\pi g^2{\vec{Q}}_1\cdot {\vec{Q}}_2+8\pi {\left(g^{\prime }\right)}^2{\vec{Q}}_1\cdot {\vec{Q}}_2+16\pi g{g}^{\prime }{\vec{Q}}_1\times {\vec{Q}}_2\right)I_3\log R_{12}+$$

$$+y^4\left(4\pi g\left(g+g^{\prime }\right)\left({\vec{Q}}_1\times {\vec{Q}}_1+{\vec{Q}}_2\times {\vec{Q}}_2\right)I_1+8{\left(g^{\prime }\right)}^2{I}_1\right)\log ϵ,$$

$$\frac{\partial g}{\partial \ell }=-16\pi (g^2+g^{\prime 2})y^4$$

$$\frac{\partial g^{\prime }}{\partial \ell }=-2\pi g{g}^{\prime }{y}^4$$

$$\frac{\partial y}{\partial \ell }=2\pi (1-g-g^{\prime })y.$$

$$\mathcal{Z}=\underset{n\to 0}{\lim }\mathrm{Tr}\exp \left(-\sum _{\mu =1}^n{\mathcal{H}}_{\mathrm{vort}}\left(Q^{(\mu )}\right)\right),$$

$${\mathcal{H}}_{\mathrm{vort}}=\sum _{ij}\sum _{\alpha \beta }{Q}_{i\alpha }{J}_{ij}^{\alpha \beta }{Q}_{j\beta }$$

$${\bar{\mathcal{Z}}}^n=\int \mathcal{D}\left[Q_{i\alpha }^{(\mu )}\right]\int \mathcal{D}\left[J_{ij}^{\alpha \beta }\right]\exp \left[-\frac{1}{2}\sum _{i,j=1}^N\sum _{\alpha ,\beta }\left(J_{ij}^{\alpha \beta }-J_0^{\alpha \beta }\right){\sigma }_{\alpha \beta }^{-2}\left(J_{ij}^{\alpha \beta }-J_0^{\alpha \beta }\right)-\sum _{\mu =1}^n\sum _{i,j=1}^N\sum _{\alpha ,\beta }\beta J_{ij}^{\alpha \beta }{Q}_{i\alpha }^{(\mu )}{Q}_{j\beta }^{(\mu )}\right].$$

$${\bar{\mathcal{Z}}}^n=const.\times \int \mathcal{D}\left[Q_{i\alpha }^{(\mu )}\right]\exp \left[\frac{1}{2}{\beta }^2\sum _{\mu ,\nu =1}^n\sum _{i,j=1}^N\sum _{\alpha ,\beta }{Q}_{i\alpha }^{(\mu )}{Q}_{i\beta }^{(\nu )}{\left({\hat{\sigma }}^2\right)}_{\alpha \beta }{Q}_{j\alpha }^{(\mu )}{Q}_{j\beta }^{(\nu )}-\beta \sum _{\mu =1}^n\sum _{i,j=1}^N\sum _{\alpha ,\beta }{J}_0^{\alpha \beta }{Q}_{i\alpha }^{(\mu )}{Q}_{j\beta }^{(\mu )}\right].$$

$$p_{\alpha }^{(\mu )}=\frac{1}{N}\sum _{i=1}^N{Q}_{i\alpha }^{(\mu )},q_{\alpha \beta }^{(\mu \nu )}=\frac{1}{N}\sum _{i,j=1}^N{Q}_{i\alpha }^{(\mu )}{Q}_{j\beta }^{(\nu )},$$

$$M(\mathbf{r})=\sum _{\alpha =1,2}{\mathbf{M}}_{\alpha }(\mathbf{r})\cos (\mathbf{n}\cdot \mathbf{r}),$$

$${\mathcal{H}}_{\mathrm{af}}=\frac{1}{2g_M}\left[{\left(\frac{1}{c_M}{\partial }_{\tau }{\mathbf{M}}_{\alpha }\right)}^2+|\nabla {\mathbf{M}}_{\alpha }{|}^2+\frac{r}{2}|{\mathbf{M}}_{\alpha }{|}^2\right]+\frac{u_0}{2}|{\mathbf{M}}_{\alpha }{|}^4-v_0{\left(|{\mathbf{M}}_1{|}^2+|{\mathbf{M}}_2{|}^2\right)}^2.$$

$${\mathcal{H}}_{\mathrm{af}}=\frac{1}{2}{\left({\partial }_{\tau }M\right)}^2+\frac{1}{2}|\nabla M{|}^2+\frac{r}{2}|M{|}^2+\frac{u}{2}{\left(|M{|}^2\right)}^2-v|{\mathbf{M}}_1{|}^2|{\mathbf{M}}_2{|}^2,$$

$${\mathbf{M}}_{\alpha }(r,\varphi )=\int d{\varphi }^{\prime }{e}^{\frac{ı}{2}\left({\varphi }^{\prime }-\varphi \right){\hat{\ell }}_3}{\mathbf{m}}_{\alpha },$$

$${\mathbf{M}}_{\alpha }=\left(\begin{array}{ccc}\cos \varphi & \mp \sin \varphi & 0\\ \pm \sin \varphi & \cos \varphi & 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{c}{m}_{1\alpha }\\ m_{2\alpha }\\ m_{3\alpha }\end{array}\right),$$

$$E_1=2\pi (|{\mathbf{m}}_X\times {\mathbf{e}}_Z{|}^2+|{\mathbf{m}}_Y\times {\mathbf{e}}_Z{|}^2)\log \Lambda =2\pi |{\mathbf{m}}_{\perp \alpha }{|}^2\log \Lambda ,$$

$$E_2=2\pi {\vec{Q}}_i\cdot {\vec{Q}}_j\left(|{\mathbf{m}}_1\times {\mathbf{e}}_Z{|}^2+|{\mathbf{m}}_2\times {\mathbf{e}}_Z{|}^2\right)\log r_{ij}=2\pi |{\mathbf{m}}_{\perp \alpha }{|}^2{\vec{Q}}_i\cdot {\vec{Q}}_j\log r_{ij}.$$

$${\mathcal{H}}_{\mathrm{vort}}=\sum _{i<j}\left(g{\vec{Q}}_i\cdot {\vec{Q}}_j+g^{\prime }{\vec{Q}}_i\times {\vec{Q}}_i\right)\log r_{ij}+\sum _i\vec{\mu }\cdot {\vec{Q}}_i.$$

$$g=m_{\perp }^2+\frac{4r+6u{m}_{\perp }^2}{\left(2v+\frac{3}{2}u{m}_{\perp }^2+\frac{v}{2}{m}_{\perp }^2\right)\left(2r+\frac{3}{2}u{m}_{\perp }^2-\frac{v}{2}{m}_{\perp }^2\right)}$$

$$g^{\prime }=-\frac{4v{m}_{\perp }^2}{\left(2v+\frac{3}{2}u{m}_{\perp }^2+\frac{v}{2}{m}_{\perp }^2\right)\left(2r+\frac{3}{2}u{m}_{\perp }^2-\frac{v}{2}{m}_{\perp }^2\right)}$$

$${\mu }_{\alpha }=\frac{1}{2}{m}_{\perp }{m}_z,$$

$$\frac{\delta \mathcal{Z}}{\mathcal{Z}}$$

Full text

I Introduction

Nonlinear and pattern-forming systems patt1 ; pattbook ; korn have numerous analogies with strongly correlated systems encountered in condensed matter physics ft1 ; ft2 , and on the methodological level they are both united through the language of field theory, which has become the standard language to describe strongly correlated electrons fradkin ; tsvelik as well as nonlinear dynamical systems cvetbook . In the field of pattern formation, some connections to condensed matter systems have been observed, see e.g. ft1 . More recently, extensive field-theoretical studies of laser systems were performed, e. g. conti1 ; conti2 ; conti3 ; conti4 , and also compared to experiment ghofraniha . However, this topic is far from exhausted and we feel many analogies between quantum many-body systems and pattern-formation dynamics remain unexplored and unexploited. In particular, nonlinear optical systems and photonic lattices are flexible and relatively cheap to build korn and they can be used to ”simulate” a broad spectrum of phenomena concerning band structure, spin ordering and conduction in strongly correlated electron systems; some of the work in this direction can be found in primeri1 ; primeri2 .

Our goal is to broaden the connections between the strongly correlated systems and nonlinear optics and to put to work the mighty apparatus of field theory to study the patterns in a nonlinear optical system from the viewpoint of phase transition theory: pattern dynamics in certain cases shows critical behavior which is analogous to phenomena seen in magnetic systems. To that end, we use the formalism of perturbative field theory and renormalization group analysis but we also perform numerical simulations from the first principles, i.e. directly integrating the equations of motion to provide an independent check of our main conclusions. We also establish a connection to an $O(3)$ antiferromagnetic model which is encountered in the study of strongly correlated electron systems. The analogy is not just qualitative: we construct the phase diagrams of both systems and find they have the same structure. Introducing disorder into the system further enriches the physics, and it is physically motivated: in optics, disorder is rooted in the imperfections of the photonic lattice, and in magnetic systems it comes from the quenched spin impurities which are regularly found in realistic samples. It turns out that in both cases a glassy phase arises. This is another important research topic and it is again appealing to realize glasses in photonic lattice systems, where the parameters are easy to tune.

I.1 On topology and vortices

The key phenomenon which governs the phenomenology of the systems studied is the existence of topologically nontrivial solutions or topological solitons raja : these are the solutions which map the physical boundary of the system to the whole configuration space of the field, so one explores all field configurations by ”going around the system”. For example, in a two-dimensional system (in the $x-y$ plane) with $U(1)$ phase symmetry, the configuration space is a circle (the phase lies between [math] and $2\pi$) and the boundary of the physical space (i.e. the two-dimensional plane) is again a circle, the ”boundary” of the plane at infinity. The topological soliton is a pattern of the $U(1)$ field which spans the whole phase circle (its phase goes from [math] to $2\pi$), as one moves around the far-away circle in the $x-y$ plane. Of course, this is the vortex – the most famous and best studied topological configuration. Similar logic leads to the classification of topological defects of other, more complicated symmetry groups. A potential source of confusion is that in nonlinear dynamics and theory of partial differential equations, the ”integrable” solutions, i.e. linearly (often also nonlinearly) stable solutions which can be obtained by inverse scattering or similar methods and which propagate through each other without interacting, are also called solitons, or more precisely dynamical solitons. In optics, they are often called spatial solitons. Dynamical solitons in nonlinear optics are a celebrated and well-studied topic rev ; kivshar1 ; kivshar2 ; kivshar3 ; newsol1 ; newsol2 , they show an amazing variety of patterns and phenomena like localization, Floquet states primeri1 , etc. But in general they do not have a topological charge. In contrast, topological solitons carry a topological charge (winding number for vortices) and their stability is rooted in topological protection (conservation of topological charge).

The phenomenon of vortices is perhaps best known in three spatial dimensions. The phase of the wavefunction can wind, forming a vortex line. These vortices are stable when the phase symmetry is broken by magnetic field. Famously, vortices may coexist with the superconducting order ($U(1)$ symmetry breaking) in type-II superconductors or exist only in the normal phase, upon destroying the superconductivity (type I). The primary example in two spatial dimensions is the vortex unbinding (BKT) phase transition of infinite order found by Berezinsky, Kosterlitz and Thouless for the planar XY model kosterthoul . The formal difference between the two- and three-dimensional vortices is that the latter give rise to an emergent gauge field; this does not happen in the XY-like system in two dimensions kleinert . While the nonlinear optical system we study is three-dimensional, its geometry and relaxational dynamics make it natural to treat it as a $2+1$-dimensional system (the $x$ and $y$ coordinates are spatial dimensions, the $z$-direction has the formal role of time, and physical time $t$ has the role of a parameter). We therefore have a similar situation to the XY model: point-like vortices in the plane (and no gauge field).

Vortex matter is known to emerge in liquid helium liqhel , Bose-Einstein condensates bec and magnetic systems vortcupsach . The basic mechanisms of vortex dynamics are thus well known. However, novel physics can arise if the system has multiple components and each of them can form vortices which mutually interact. This is precisely our situation – we have a system of two laser beams propagating in opposite directions, and we will compare it to a two-component antiferromagnet. So far, such situations have been explored in multi-component superconductors multisc which have attracted some attention, as they can be realized in magnesium diboride multiscmg . But these are again bulk systems, not planar. Vortices in planar multi-component systems have not been very popular, an important exception being the two-component Bose-Einstein condensates of becmulti , which were found to exhibit complex vortex dynamics; in these systems, contrary to our case, the two components have an explicit attractive interaction, unlike our case where they interact indirectly, by coupling to the total light intensity (of both components).

I.2 The object of our study

In this paper we study phases and critical behavior of topological configurations (vortices and vortex lattices) in a specific and experimentally realizable nonlinear optical system: laser beams counterpropagating (CP) through a photorefractive (PR) crystal. This means we have an elongated PR crystal (with one longitudinal and two transverse dimensions) and two laser beams shone onto each end. We thus effectively have two fields, one forward-propagating and one backward propagating. The optical response of the crystal depends nonlinearly on the total intensity of both beams, which means the beams effectively interact with each other. This system has been thoroughly investigated for phenomena such as dynamical solitons prlpet ; stabanal ; rev , vortex stability on the photonic lattice vortlat2pet ; prapet ; vortprepet ; kivshar1 ; kivshar2 ; kivshar3 ; kivshar4 and global rotation gaussrotpet . We will see that the CP beams are an analog of the two-component planar antiferromagnet, which can further be related to some realistic strongly-correlated materials sachbook ; jurms0 ; jurms2 . The two beams are now equivalent to two sublattices which interact through a lattice deformation or external field. The PR crystal is elongated and the axial propagation direction has the formal role of time, which has a finite span, the length of the crystal. For the antiferromagnet, the third axis is the usual imaginary time compactified to the radius $1/T$, i.e. inverse temperature. Both systems contain vortices as topological defects, i.e. solutions with integer topological charge. In the PR optical system, vortices arise as a consequence of the $U(1)$ symmetry of the electromagnetic field. In antiferromagnets we consider, the $O(3)$ symmetry of the antiferromagnet gives rise to $\mathbb{Z}_{2}$-charged defects, which exhibit the same interactions as the vortices. The optical system is not subject to noise, i.e. it lives at zero temperature, thus the criticality we talk about is obviously not the same as thermodynamic phase transitions. Phase transitions happen upon varying the parameters, not temperature, so they may be described as quantum critical phenomena in the broad sense taken in sachbook – any critical behavior controlled not by thermal fluctuations but by parameter dependence.

In the PR counterpropagating beam system, our focus are the vortices but in order to study them we need to do some preparational work. We first recast the system in Lagrangian and then in Hamiltonian form so it can be studied as a field theory, which depends parametrically on the time $t$. Then we consider the time dynamics of the system and show that in a broad parameter range the patterns relax to a static configuration which can be studied within equilibrium field theory. Along the way, we also study the stability of topologically trivial (vortex-free) configurations and then consider the phases of the static vortex configurations. The analytical insight we obtain also allows us to avoid overextensive numerics – analytical construction of the phase diagram tells us which patterns can in principle be expected in different corners of the parameter space. By ”blind” numerical approach this result could only be found through many runs of the numerics.

In the antiferromagnetic spin system the non-topological excitations are simple – they are spin waves, perturbed away from the noninteracting solution by the quartic terms in the potential. There are no dynamical solitons. But we will see that topological excitations lead to a phase diagram which, after reasonable approximations, can be exactly mapped to the phase diagram of the photorefractive crystal. The reason is that both can be reduced to an effective Hamiltonian for a two-component vortex system, i.e., every vortex has two charges, or two ”flavors”. In the photorefractive crystal it happens naturally, as there are two beams, forward- and backward propagating. In the Heisenberg antiferromagnet it is less obvious, and is a crucial consequence of the collinearity of the spin pattern. We will focus on common properties of the two systems and map the phase diagrams onto each other. In the antiferromagnetic system, different phases are separated by quantum phase transitions – phase transitions driven by the quantum fluctuations instead of temperature.

I.2.1 On disorder

It is known that impurities pin the vortices and stabilize them. This leads to frozen dynamics even though no symmetry is broken, the phenomenon usually associated with glasses. In simple systems such as the Ising model with disorder one generically has two phases: the disordered (paramagnetic) phase remains and the ordered (magnetic) phase is replaced by a regime with algebraic correlations and no true order. In many cases, such phases are called glasses. The exact definition of a glass is lacking; normally, they show (i) long-range correlations (ii) absence of long-range order, i.e. of a nonzero macroscopic order parameter (iii) ”frozen dynamics”, i.e. free energy landscape with numerous local minima in which the system can spend a long time parisi ; spinglass . While the most popular example are probably spin glasses in Ising-like models such as Sherington-Kirkpatrick and Edwards-Anderson model, glasses are also known to appear in the XY model with disorder in two dimensions, the Cardy-Ostlund model, which postulates both random couplings and a random magnetic field xyzoran ; xy1 ; xy2 . Our model is essentially a two-flavor generalization of the XY model, although in order to solve it we need to simplify it. According to xyzoran ; xy1 ; xy2 the details differ depending on how the disorder is implemented but the two-phase system (paramagnetic i.e. disordered, and glass) is ubiquitous. In the two-component version, the phase diagram becomes richer, and on top of the glassy phase and the insulator (disordered) phase we find a few other phases. In nonlinear optics, the topic of random lasers has attracted considerable attention conti1 ; conti2 ; conti3 ; conti4 ; conti5 . Here one has a complex version of the XY model, with the additional complication that not only phase but also amplitude is free to vary, but only with random couplings (no random field). On top of the glassy and the disordered phase, one or two additional phases appear.

In the presence of disorder the relation to magnetic systems in condensed matter physics is very inspiring, since a number of complex materials show different ordering mechanisms (spin and charge density waves, superconductivity, etc) in parallel with significant influence of disorder. Just as in the disorder-free case, we are particularly interested in possible spin-glass phenomena in doped insulating $O(3)$ antiferromagnets jurms0 ; jurms1 ; jurms2 ; jurms3 ; jurms4 and in the last section we will discuss also the spin-glass phase in such systems.

I.3 The plan of the paper

The structure of the paper is as follows. In the next section we describe the dynamical system which lies at the core of this paper: counterpropagating laser beams in a photorefractive crystal. We give the equations of motion and repackage them in the Lagrangian form. In section III we study the vortex dynamics: we construct the vortex Hamiltonian and classify the order parameters. Then we study the renormalization group (RG) flow and obtain the phase diagram. Finally, we discuss the important question of how to recognize the various phases in experiment: what do the light intensity patterns look like and how they depend on the tunable parameters. Section IV brings the same study for the system with disorder. After describing the disordered system, we perform the replica trick for the disordered vortex Hamiltonian and solve the saddle-point equations to identify the phases and order parameters, again refining the results with RG calculations. The fifth section takes a look at a doped collinear antiferromagnet, a model encountered in the description of many strongly coupled materials, and shows how the dynamics of topological solitons is again described by a two-flavor vortex Hamiltonian. We discuss the relation between the phase diagrams of the two systems and the possibilities of modeling the condensed matter systems experimentally by the means of photorefractive optics. The last section sums up the conclusions. In Appendix A we describe the numerical algorithm we use to check the analytical results for the phase diagram. In Appendix B we show in detail that the CP beams are capable of reaching equilibrium (i.e. stop changing in time)– if they would not, the application of equilibrium field theory would not be justified. Appendix C discusses the stability of nonvortex configurations – although somewhat peripheral to the main topic of the paper, it is useful to better understand the geometry of patterns. In Appendix D we give the (routine) algebra that yields the vortex interaction Hamiltonian from the microscopic equations. Appendix E contains an improved mean-field theory for the clean system, which we do not use much throughout the paper but we include it for completeness (we prefer either the simplest single-vortex mean-field reasoning, or the full RG analysis, which are described in the main text). Appendix F discusses an important technicality concerning the CP geometry, i.e. the specific boundary conditions of the CP beam system where the boundary conditions for one beam are given at the front face and for the other at the back face of the crystal. Appendix G contains some details on mean-field and RG calculations of the phase diagram for the dirty system: the dirty case includes some tedious algebra we feel appropriate to leave out from the main text.

II The model of counterpropagating beams in the photorefractive crystal

We consider a photorefractive crystal of length $L$ irradiated by two laser beams. The beams are paraxial and propagate head-on from the opposite faces of the crystal in the $z$-direction. Photorefractive crystals induce self-focusing of the beams – the vacuum (linear) wave equation is modified by the addition of a friction-like term, so the diffusion of the light intensity (the broadening of the beam) is balanced out by the convergence of the beam onto an ”attractor region”. The net result is the balance between the dissipative and scattering effects, allowing for stable patterns to form. The physical ground for this is the redistribution of the charges in the crystal due to the Kerr effect. The nonlinearity, i.e. the response of the crystal to the laser light is contained in the change of the refraction index which is determined by the local charge density. A sketch of the system is given in Fig. 1. Before entering the crystal, the laser beams can be given any desirable pattern of both intensity and phase. In particular, one can create vortices (winding of the phase) making use of the phase masks korn or other, more modern ways.

Assuming the electromagnetic field of the form $\mathbf{E}=e^{\imath\omega t+\imath\mathbf{q}\cdot\mathbf{r}}\left(Fe^{ikz}+Be^{-ikz}\right)$, we can write equations for the so-called envelopes $F$ and $B$ of the forward- and backward propagating beams along the $z$-axis (the frequency, transverse and longitudinal momentum are denoted respectively by $\omega,\mathbf{q},k$). The wave equations for $F$ and $B$ are now:

[TABLE]

where the plus and minus signs on the left-hand side stand for the forward- and backward-propagating component of the beam amplitude doublet $\Psi\equiv(\Psi_{+},\Psi_{-})\equiv(F,B)$, and $\Gamma$ is the dimensionless PR coupling constant. The two beams (flavors of the field $\Psi$) will from now on be denoted either by $F/B$ or more often by $\Psi_{\pm}$. We will use $\alpha$ as the general flavor index for summation, e.g. $\Psi_{1\alpha}\Psi_{2\alpha}=\Psi_{1+}\Psi_{2+}+\Psi_{1-}\Psi_{2-}$. The charge field $E$ on the right-hand side of the equation is the electric field sourced by the charges in the crystal (i.e., it does not include the external electric field of the beams). Its evolution is well represented by a relaxation-type equation rev :

[TABLE]

Here, $I\equiv I_{\Psi}+I_{x}$ is the total light intensity at a given point, $I_{\Psi}\equiv|F|^{2}+|B|^{2}$ is the beam intensity and $I_{x}$ the intensity of the fixed background. The meaning of $I_{x}$ is that the crystal is all the time irradiated by some constant light source, independent of the counter-propagating beams with envelopes $F,B$. We will usually take a periodic lattice as the background, allowing also for the defects (missing cells) in the lattice when studying the effects of disorder. The relaxation time is $\tau$. The time derivative $\partial_{t}E$ is divided by $1+I$, meaning that the polarizability of the crystal depends on the total light intensity: strongly irradiated regions react faster. In the numerical calculations, we solve Eqs. (1-2) with no further assumptions, as explained in Appendix A. For analytical results we will need to transform them further assuming a vortex pattern.

The equation for the charge field has no microscopic basis, it is completely phenomenological, but it excellently represents the experimental results korn . Notice that the derivative $\partial_{t}E$ in (2) is strictly negative (since intensity is non-negative): it thus has the form of a relaxation equation, and one expects that a class of solutions exists where $\partial_{t}E(t\to\infty)\to 0$, i.e. the system relaxes to a time-independent configuration. We show this in Appendix B; in the main text we will not discuss this issue but will simply take the findings of the Appendix B for granted. Notice that there are also parameter values for which no equilibrium is reached fisher ; fisher2 ; gaussrotpet .

For slow time evolution (in the absence of pulses), we can Laplace-transform the equation (2) in time ($E(t)\mapsto E(u)=\int_{0}^{\infty}dte^{-ut}E(t)$) to get the algebraic relation

[TABLE]

The original system (1) can now be described by the Lagrangian:

[TABLE]

where $\sigma_{3}$ is the Pauli matrix $\sigma_{3}=\mathrm{diag}(1,-1)$. One can introduce the effective potential

[TABLE]

so we can write the Lagrangian as $\mathcal{L}=\imath\Psi^{\dagger}\sigma_{3}\partial_{z}\Psi-|\nabla\Psi|^{2}-V_{eff}(\Psi^{\dagger},\Psi)$. This is the Lagrangian of a non-relativistic field theory (a nonlinear Schrödinger field equation) in $2+1$ dimensions $(x,y;z)$, where the role of time is played by the longitudinal distance $z$, and the physical time $t$ (or $u$ upon the Laplace transform) is a parameter. The span of the $z$ coordinate $0<z<L$ will influence the behavior of the system, while the dimensions of the transverse plane are not important for the effects we consider.

Our main story is now the nature and interactions of the topologically nontrivial excitations in the system (4). A task which is in a sense more basic, the analysis of the topologically trivial vacua of (4) and perturbative calculation of their stability, is not of our primary interest now, in part because this was largely accomplished by other methods in stabanal ; prlpet . We nevertheless give a quick account in Appendix C; first, because some conclusions about the geometry of the patterns can be carried over to vortices, and second, to give another example of applying the field-theoretical formalism whose power we wish to demonstrate and popularize in this paper.

III Vortices and mean field theory of vortex interactions

III.1 The classification of topological solutions and the vortex Hamiltonian

Now we discuss the possible topological solitons in our system. Remember once again that they differ from dynamical solitons such as those studied in rev and references therein. In order to classify the topologically nontrivial solutions, consider first the symmetries of the Lagrangian (4). It describes a doublet of 2D complex fields which interact solely through the phase-invariant total intensity $I=\Psi^{\dagger}\Psi$ (and the spatial derivative term $|\nabla\Psi|^{2}$), while in the kinetic term $\Psi^{\dagger}\sigma_{3}\partial_{z}\Psi$ the two components have opposite signs of the ”time” derivative, so this term cannot be reduced to a functional of $I$. The intensity $I$ has the symmetry group $SU(2)$ (the isometry group of the three-dimensional sphere in Euclidean space) and the kinetic term has the group $SU(1,1)$ (the transformations which leave the combination $|F|^{2}-|B|^{2}$ invariant, i.e. the isometry of the hyperboloid). The intersection of these two is the product $U(1)_{F}\otimes U(1)_{B}$: the forward- and backward-propagating doublet $(F,B)$ has phases $\theta_{F,B}$ which can be transformed independently, as $\theta_{F,B}\mapsto\theta_{F,B}+\delta\theta_{F,B}$.

The classification of possible topological solitons is straightforward from the above discussion mermin . They can be characterized in terms of homotopy groups. To remind, the homotopy group $\pi_{n}$ of the group $G$ is the group of transformations which map the group manifold of $G$ onto the $n$-dimensional sphere $\mathbb{S}_{n}$. In $D$-dimensional space the group $\pi_{D-1}$ therefore classifies what a field configuration looks like from far away (from infinity): it classifies the mappings from the manifold of the internal symmetry group of the system to the spherical ”boundary shell” in physical space at infinity. Since the beams in our PR crystal effectively see a two-dimensional space (we regard $z$ as time), we need the first homotopy group $\pi_{1}$ to classify the topological solitons. Since $\pi_{1}(U(1))=\pi_{1}(\mathbb{S}_{1})=\mathbb{Z}$ and $\pi_{1}(\mathbb{G}\otimes\mathbb{G})=\pi_{1}(\mathbb{G})\otimes\pi_{1}(\mathbb{G})$ for any group $\mathbb{G}$, the topological solutions are flavored vortices, and the topological charge is the pair of integers $\{Q_{F},Q_{B}\}$.

Let us now derive the effective interaction Hamiltonian for the vortices and study the phase diagram. In principle, this story is well known: for a vortex at $\mathbf{r}_{0}$, in the polar coordinates $(r,\phi)$, we write $\Psi(\mathbf{r})=\psi\exp\left(\imath\theta\left(\mathbf{r}\right)\right)$ for $|\mathbf{r}-\mathbf{r}_{0}|/|\mathbf{r}_{0}|\ll 1$, and a vortex of charge $Q$ has $\theta(\phi)=Q\phi/2\pi$. In general the phase has a regular and a singular part, $\nabla\Psi=\psi(\nabla\delta\theta+\nabla\times\zeta\mathbf{e}_{z})$, where finally $\zeta=Q\log|\mathbf{r}-\mathbf{r}_{0}|$. The difference in the CP beam system lies in the existence of two beam fields (flavors) and the non-constant amplitude field $\psi_{\pm}(r)$, so the vortex looks like

[TABLE]

Inserting this solution into the equations of motion (or, equivalently, the Lagrangian) it is just a matter of algebra to obtain the vortex Hamiltonian, analogous to the well-known one but with two components (flavors) and their interaction. We refer the reader to the Appendix D for the full derivation. The outcome is perhaps expected: we get the straightforward generalization of the familiar Coulomb gas picture for the XY model where all interactions of different flavors – $F$-$F$, $B$-$B$ and $F$-$B$ are allowed. In order to write the Hamiltonian (and further manipulations with it) in a concise way, it is handy to introduce shorthand notation $\vec{Q}\equiv(Q_{+},Q_{-})$, $\vec{Q}_{1}\cdot\vec{Q_{2}}\equiv Q_{1+}Q_{2+}+Q_{1-}Q_{2-}$, and $\vec{Q}_{1}\times\vec{Q}_{2}\equiv Q_{1+}Q_{2-}+Q_{1-}Q_{2+}$. For the self-interaction within a vortex $\vec{Q}_{1}$, we have $\vec{Q}_{1}\cdot\vec{Q}_{1}=Q_{1+}^{2}+Q_{1-}^{2}$ but $\vec{Q}_{1}\times\vec{Q}_{1}\equiv Q_{1+}Q_{1-}$ (i.e. there is a factor of $2$ mismatch with the case of two different vortices). Now for vortices at locations $\mathbf{r}_{i},i=1,\ldots,N$ with charges $\{Q_{i+},Q_{i-}\}$ we get:

[TABLE]

The meaning of the Hamiltonian (7) is obvious. The first term is the Coulomb interaction of vortices; notice that only like-flavored charges interact through this term (because the kinetic term $|\nabla\Psi|^{2}$ is homogenous quadratic). The second term is the forward-backward interaction, also with Coulomb-like (logarithmic) radial dependence. This interaction comes from the mixing of the $F$- and $B$-modes in the fourth term in Eq. (53), and it is generated, as we commented in Appendix D, when the amplitude fluctuations $\delta\psi_{\alpha}(r)$, which couple linearly to the phase fluctuations, are integrated out. In a system without amplitude fluctuations, i.e. classical spin system, this term would not be generated. The third and fourth terms constitute the energy of the vortex core. The self-interaction constants $g_{0},g_{1}$ are of course dependent on the vortex core size and behave roughly as $g\log a/\epsilon,g^{\prime}\log a/\epsilon$, where $\epsilon$ is the UV cutoff. The final results will not depend on $\epsilon$, as expected, since $g_{0},g_{1}$ can be absorbed in the fugacity $y$ (see the next subsection). Expressions for the coupling constants in terms of original parameters are given in (62).

In three space dimensions, vortices necessitate the introduction of a gauge field kleinert which, in multi-component systems, also acquires the additional flavor index multisc ; multi15 . In our case, there is no emergent gauge field and the whole calculation is a rather basic exercise at the textbook level but the results are still interesting in the context of nonlinear optics and analogies to magnetic systems: they imply that the phase structure (vortex dynamics) can be spotted by looking at the intensity patterns (light intensity $I$ or local magnetization $\mathcal{M}$, see the penultimate section).

III.2 The phase diagram

III.2.1 The mean-field theory for vortices

The phases of the system can be classified at the mean field level, following e.g. parisi ; kleinert . In order to do that, one should construct the partition function, assuming that well-defined time-independent configuration space exists. We have already mentioned the question of equilibration, and address it in detail in Appendix B. Knowing that the system reaches equilibrium (in some part of the parameter space), we can count the ways in which a system of vortices can be placed in the crystal – this is by definition the partition function $\mathcal{Z}$. First, the number of vortices $N$ can be anything from [math] to infinity, second, the vortex charges can be arbitrary, and finally the number of ways to place each vortex in the crystal is simply the total surface section of the crystal divided by the size of the vortex. Then, each vortex carries a Gibbs weight proportional to the energy, i.e. the vortex Hamiltonian (7) for a single vortex.111Again, this is not generally true for out-of-equilibrium configurations but if the system reaches equilibrium, i.e., stable fixed point, this follows by usual statistical mechanics reasoning. Let us focus first on a single vortex. If the vortex core has linear dimension $a$ and the crystal cross section linear dimension $\Lambda$, the vortex can be placed in any of the $(\Lambda/a)^{2}$ cells (and in the mean-field approach we suppose the vortex survives all the way along the crystal, from $z=0$ to $z=L$, so there is no additional freedom of placing it along some subinterval of $z$). This gives

[TABLE]

Remember that $\mathcal{H}$ is energy density along the $z$ axis, so it appears multiplied by $L$. The factor $\log(\Lambda/a)$ in the second term of the exponent comes from the Coulomb potential of a single vortex (in a plane of size $\Lambda$). The exponent can be written as $-L\mathcal{F}^{(1)}$, with $\mathcal{F}^{(1)}=\mathcal{H}_{1}-(1/L)S_{1}$, recovering the relation between the free energy $\mathcal{F}^{(1)}$ and entropy $S^{(1)}$ of a single vortex. The entropy comes from the number of ways to place a vortex of core size $a$ in the plane of size $\Lambda\gg a$: $S\sim\log(\Lambda/a)^{2}$. Suppose for now that elementary excitations have $|Q_{\pm}|\leq 1$, as higher values increase the energy but not the entropy, so they are unlikely (when only a single vortex is present). Now we can consider the case of single-charge vortices with possible charges $(1,0),(-1,0),(0,1),(0,-1)$, and the case of two-charge vortices where $F$- and $B$-charge may be of the same sign or of the opposite sign, $(1,1),(-1,-1),(1,-1),(-1,1)$:

[TABLE]

Now we identify four regimes, assuming that $g,g^{\prime}>0$:222One specificity of multi-component vortices is that the coupling constants may be negative, as can be seen from (62). In that case the ordering of the four regimes (how they follow each other upon dialing $L$) changes but the overall structure remains.

For $L>2/g$ a vortex always has positive free energy so vortices are unstable like in the low-temperature phase of the textbook BKT system. This is the vortex-free phase where the phase $U(1)_{F}\otimes U(1)_{B}$ does not wind. This phase we logically call vortex insulator in analogy with the single-flavor case. 2. 2.

For $2/g>L>1/(g-g^{\prime}/2)$ a double-flavor vortex always has positive free energy but single-flavor vortices are stable; in other words, there is proliferation of vortices of the form $\vec{Q}=(Q_{+},0)$ or $\vec{Q}=(0,Q_{-})$. This phase is like the conductor phase in a single-component XY model, and the topological excitations exist for the reduced symmetry group, i.e. for a single $U(1)$. We thus call it vortex conductor; it is populated mainly by single-flavor vortices $(Q,0)$, $(0,Q)$. 3. 3.

For $1/(g-g^{\prime}/2)>L>1/(g+g^{\prime}/2)$ double-vortex formation is only optimal if the vortex has $Q_{+}+Q_{-}=0$ which corresponds to the topological excitations of the diagonal $U(1)_{d}$ symmetry subgroup, the reduction of the total phase symmetry to the special case $(\theta_{F},\theta_{B})\mapsto(\theta_{F}+\delta\theta,\theta_{B}-\delta\theta)$. In other words, vortices of the form $(Q_{+},-Q_{+})$ proliferate. Here higher charge vortices may be more energetically favorable than unit-charge ones, contrary to the initial simplistic assumption, the reason being that the vortex core energy proportional to $gQ_{+}^{2}$ may be more than balanced out by the intra-vortex interaction proportional to $-g^{\prime}Q_{+}^{2}$ (depending on the ratio of $g$ and $g^{\prime}$). This further means that there may be multiple ground states of equal energy (frustration). We thus call this case frustrated vortex insulator (FI); it is populated primarily with vortices of charge $(Q,-Q)$. 4. 4.

For $1/(g+g^{\prime}/2)>L$ vortex formation always reduces the free energy, no matter what the relation between $Q_{+}$ and $Q_{-}$ is, and each phase can wind separately: $(\theta_{F},\theta_{B})\mapsto(\theta_{F}+\delta\theta_{F},\theta_{B}+\delta\theta_{B})$. Vortices of both flavors proliferate freely at no energy cost and for that reason we call this phase vortex perfect conductor (PC). We deliberately avoid the term superconductor to avoid the (wrong) association of this phase with the vortex lines and type I/type II superconductors familiar from the 3D vortex systems: remember there is no emergent gauge field for the vortices in two spatial dimensions, and we only have perfect conductivity in the sense of zero resistance for transporting the (topological) charge, but no superconductivity in the sense of breaking a gauge symmetry.

A more systematic mean-field calculation will give the phase diagram also for an arbitrary number of vortices. This is not so interesting as it already does not require much less work than the RG analysis, which is more rigorous and more accurate for this problem. For completeness, we give the multi-vortex mean-field calculation in Appendix E.

One might worry that the our whole approach approach misses the CP geometry of the problem, i.e. the fact that the $\Psi_{+}$ field has a source at $z=0$ and the $\Psi_{-}$ field at $z=L$. In Appendix F we show that nothing is missed at the level of approximations taken in this paper, i.e. mean-field theory in this subsubsection and the lowest-order perturbative RG in the next one. Roughly speaking, it is because the sources are irrelevant in the RG sense – the bulk configuration dominates over the boundary terms. The Appendix states this in much more precise language.

III.2.2 RG analysis

We have classified the symmetries and thus the phases of our system at the mean-field level. To describe quantitatively the borders between the phases and the phase diagram, we will perform the renormalization group (RG) analysis. Here we follow closely the calculation for conventional vortex systems kleinert . We consider the fluctuation of the partition function $\delta\mathcal{Z}$ upon the formation of a virtual vortex pair at positions $\mathbf{r}_{1},\mathbf{r}_{2}$ with charges $\vec{q},-\vec{q}$, (with $\mathbf{r}_{1}+\mathbf{r}_{2}=2\mathbf{r}$ and $\mathbf{r}_{1}-\mathbf{r}_{2}=\mathbf{r}_{12}$), in the background of a vortex pair at positions $\mathbf{R}_{1},\mathbf{R}_{2}$ (with $\mathbf{R}_{1}+\mathbf{R}_{2}=2\mathbf{R}$ and $\mathbf{R}_{1}-\mathbf{R}_{2}=\mathbf{R}_{12}$) with charges $\vec{Q}_{1},\vec{Q}_{2}$. This is a straightforward but lengthy calculation and we state just the main steps. First, it is easy to show that the creation of single-charge vortices is irrelevant for the RG flow so we disregard it. Also, we can replace the core self-interaction constants $g_{0,1}$ with the fugacity parameter defined as $y\equiv\exp\left[-\beta\left(g_{0}+g_{1}\right)\log\epsilon\right]$. Here we introduce the notation $\beta\equiv L$ in analogy with the inverse temperature $\beta$ in standard statistical mechanics, in order to facilitate the comparison with the literature on vortices in spin systems, and also with antiferromagnetic systems in Section V.333Of course, the physical meaning of $\beta$ in our system is very different: we have no thermodynamic temperature or thermal noise, and the third law of thermodynamics is not satisfied for the ”temperature” $1/\beta=1/L$. We merely use the $\beta$-notation to emphasize the similarity between free energies of different systems, not as a complete physical analogy.

Now from the vortex Hamiltonian $\mathcal{H}_{\mathrm{vort}}$ the fluctuation equals (at the quadratic order in $y$ and $r$):

[TABLE]

Notice that $\nabla$ is taken with respect to $\mathbf{r}$. The above result is obtained by expanding the Coulomb potential in $r_{12}$ (the separation between the virtual vortices being small because of their mutual interaction) and then expanding the whole partition function (i.e. the exponent in it) in $y$ around the equilibrium value $\mathcal{Z}$. The term depending on the separation $r_{12}$ is the mutual interaction energy of the virtual charges, and the subsequent term proportional to $r^{2}$ is the interaction of the virtual vortices with the external ones (the term linear in $r$ cancels out due to isotropy). Then by partial integration and summation over $q_{\pm}\in\{1,-1\}$ we find

[TABLE]

with $I_{n}=\int_{\epsilon a}^{\Lambda a}drr^{n+g+g^{\prime}}$. Now, taking into account the definition of the fugacity $y$, rescaling $\Lambda\mapsto\Lambda(1+\ell)$, performing the spatial integrals, and expanding over $\ell$, we can equate the bare quantities $g,g^{\prime},y$ in (7) with their corrected values in $\mathcal{Z}+\delta\mathcal{Z}$ to obtain the RG flow equations:

[TABLE]

Now let us consider the fixed points of the flow equations. If one puts $g^{\prime}=0$, they look very much like the textbook XY model RG flow, except that the fugacity enters as $y^{4}$ instead of $y^{2}$ (simply because every vortex contributes two charges). They yield the same phases as the mean-field approach as it has to be, but now we can numerically integrate the flow equations to find exact phase borders. The fugacity $y$ can flow to zero (meaning that the vortex creation is suppressed and the vortices tend to bind) or to infinity, meaning that vortices can exist at finite density. At $y=0$ there is a fixed line $g+g^{\prime}=1$. This line is attracting for the half-plane $g+g^{\prime}>1$; otherwise, it is repelling. There are three more attraction regions when $g+g^{\prime}<1$. First, there is the point $y\to\infty,g=g^{\prime}=0$ which has no analog in single-component vortex systems. Then, there are two regions when $g\to\infty$ and $g^{\prime}\to\pm\infty$ (and again $y\to\infty$). Of course, the large $g,g^{\prime}$ regime is strongly interacting and the perturbation theory eventually breaks down, so in reality the coupling constants grow to some finite values $g_{*},g^{\prime}_{*}$ and $g_{**},g^{\prime}_{**}$ rather than to infinities. The situation is now the following:

The attraction region of the fixed line is the vortex insulator phase: the creation rate of the vortices is suppressed to zero. 2. 2.

The zero-coupling fixed point attracts the trajectories in the vortex perfect conductor phase: only the fugacity controls the vortices and arbitrary charge configurations can form. Numerical integration shows that this point also has a finite extent in the parameter space. 3. 3.

In the attraction region of the fixed point with $g_{*}<0$ and $g^{\prime}_{*}>0$ (formally they flow to $-\infty$ and $+\infty$, respectively), same-sign $F$- and $B$-charges attract each other and those with the opposite sign which repel each other. This is the frustrated insulator. 4. 4.

The fixed point with $g_{**},g^{\prime}_{**}<0$ (formally both flow to $-\infty$) corresponds to the conductor phase.

The RG flows in the $g-g^{\prime}$ plane are given in Fig. 2. Full RG calculation is given in panel (b); for comparison, we include also the mean-field phase diagram (following from the previous subsubsection and Appendix E) in the panel (a). In the half-plane $g+g^{\prime}>1$ every point evolves toward a different, finite point $(g,g^{\prime})$ in the same half-plane. In the other half-plane we see the regions of points moving toward the origin or toward one of the two directions at infinity. The PC phase (the attraction region of the point $(0,0)$) could not be obtained from the mean field calculation (i.e. it corresponds to the single point at the origin at the mean field level).

It may be surprising that the coupling constants can be negative, with like charges repelling and opposite charges attracting each other. However, this is perfectly allowed in our system. In the usual XY model, the stiffness is proportional to the kinetic energy coefficient and thus has to be positive. Here, the coupling between the fluctuations of $F$- and $B$-beams introduces other contributions to $g,g^{\prime}$ and the resulting expressions (62) give bare values of $g,g^{\prime}$ that can be negative, and the stability analysis of the RG flow clearly shows that for nonzero $g^{\prime}$, the flow can go toward negative values even if starting from a positive value in some parameter range. If we fix $g^{\prime}=0$, the flow equations reproduce the ones from the single-component XY model, and the phase diagram is reduced to just the $g^{\prime}=0$ line. If we additionally suppose that the bare value of $g$ is non-negative, than we are on the positive $g^{\prime}=0$ semiaxis in the phase diagram – here we see only two phases, insulator (no vortices, $g\to\mathrm{const.}$) and perfect conductor ($g\to 0$). However, for $g^{\prime}$ fixed to zero (that is, with a single flavor only), the perfect conductor reduces to the usual conductor phase of the single-component XY model – in other words, we reproduce the expected behavior.

Physically, it is preferable to give the phase diagram in terms of the quantities $\Gamma,\tau,I,I_{x},L$ that appear in the initial equations of motion (1-2): the light intensities can be directly measured and controlled, whereas the relaxation time and the coupling cannot, but at least they have a clear physical interpretation. The relations between these and the effective Hamiltonian quantities $y,g,g^{\prime}$ are found upon integrating out the intensity fluctuations to obtain (7) and the explicit relations are stated in (62). Making use of these we can easily plot the phase diagram in terms of the physical quantities for comparison with experiment. However, for the qualitative understanding we want to develop here it is much more convenient to use $g,g^{\prime}$ as the phase structure is much simpler.

As an example, we plot the $\Gamma$-$g^{\prime}$ diagram in Fig. 3 (we have kept $g^{\prime}$ to keep the picture more informative; the $\Gamma$-$L$ and $\Gamma$-$I$ diagrams contain multiple disconnected regions for each phase). The non-interacting fixed point $g=g^{\prime}=0$ is now mapped to $\Gamma=0$. The tricritical point where the PC, the FI and the conductor phases meet is at $R=1$. Therefore, the rule of thumb is that low couplings $\Gamma$ produce stable vortices with conserved charges – the perfect vortex conductor. Increasing the coupling pumps the instability up, and the kind of instability (and the resulting phase) is determined by the relative strength of the photonic lattice compared to the propagating beams. Obviously, such considerations are only a rule of thumb and detailed structure of the diagram is more complex. This is one of the main motives of this study – blind numerical search for patterns without the theoretical approach adopted here would require many runs of the numerics for a good understanding of different phases.

III.3 Geometry of patterns

Now we discuss what the intensity pattern $I(\mathbf{r})$ looks like in various phases, for various boundary conditions. This is very important as this is the only thing which can be easily measured in experiment – phases $\theta_{\alpha}$ are not directly observable, while the intensity distribution is the direct outcome of the imaging of the crystal prlpet . We shall consider three situations. The first is a single Gaussian beam on zero background ($I_{x}=0$), with Gaussian initial intensity profile $|F(z=0,\mathbf{r})|^{2}=|B(z=L,\mathbf{r})|^{2}=\mathcal{N}\exp(-r^{2}/2s^{2})$ and possibly non-zero vortex charges: $\mathrm{arg}\Psi_{\pm}(\mathbf{r})\sim\exp\left(Q_{F,B}\phi\right)$, with $\mathbf{r}=(r\cos\phi,r\sin\phi)$. The second case is a quadratic vortex lattice of $F$- and $B$-beams, so the initial beam intensity is $I_{0}=\sum_{i,j}\exp\left(-\left(x-x_{i}\right)^{2}/2s_{0}^{2}-\left(y-y_{i}\right)/2s_{0}^{2}\right)$), with $x_{i+1}-x_{i}=y_{i+1}-y_{i}\equiv b=\mathrm{const.}$, the situation particularly relevant for analogies with condensed matter systems. In the third case we have again a quadratic vortex lattice but now on top of the background photonic square lattice, which is either coincident or off-phase (shifted for half a lattice spacing) with the beam lattice. The background intensity is thus of the form $I_{x}=\sum_{i,j}\exp\left(-\left(x-x_{i}\right)^{2}/2s^{2}-\left(y-y_{i}\right)/2s^{2}\right)$.

First of all, it is important to notice that there are two kinds of instabilities that can arise in a vortex beam:444They are distinct from the bifurcations which happen also in topologically trivial beam patterns and lead to the instability which eventually destroys optical (non-topological) solitons. These instabilities have been analyzed in the Appendix C and in more detail in stabanal , where the authors have found them to start from the edge of the beam and result in the classical ”walk through the dictionary of patterns”.

• •

There is an instability which originates in the disbalance between the diffusion and self-focusing (crystal response) in favor of diffusion in high-gradient regions: if a pattern $I(x,y)$ has a large gradient $\nabla I$, the kinetic term in the Lagrangian (4), i.e. the diffusion term in (1) is large and the crystal charge response is not fast enough to balance it as we travel along the $z$-axis, so the intensity rapidly dissipates and the pattern changes. Obviously, the vortex core is a high-gradient region so we expect it to be vulnerable to this kind of instability. This is indeed the case: in the center of the vortex the intensity diminishes, a dark region forms and the intensity moves toward the edges. We dub this the core or central instability (CI), and in the effective theory it can be understood as the decay of states with low fugacity $y$, i.e. high self-interaction constants $g_{0},g_{1}$. This instability prevents the formation of vortices in the insulator phase, or limits it in the frustrated insulator and conductor phases.

• •

There is an instability stemming from the dominance of diffusion over self-focusing in low-intensity regions of sufficient size and/or convenient geometry. At low intensity, the charge response is nearly proportional to $I$ (from Eq. 2), so if $I$ is small diffusion wins and the intensity dissipates. If there is sufficient inflow of intensity from more strongly illuminated regions, it may eventually balance the diffusion; but if the pattern has a long ”boundary”, i.e. outer region of low intensity, it will not happen and the pattern will dissipate out, or reshape itself to reduce the low-intensity region. We call this case the edge instability (EI). For a vortex, it happens when the positive and negative vortex charges tend to redistribute due to Coulomb attraction and repulsion. In our field theory Hamiltonian (7), this instability dominates in the conductor and perfect conductor phases.

Let us first show how the CI and EI work for a single beam with nonzero vortex charge. In Fig. 4 we show the intensity patterns for a single vortex with charges $(1,0)$ and $(3,0)$ as the $x-y$ cross sections (transverse profiles) in the middle of the crystal, i.e. for $z=L/2$. The parameters chosen ($\Gamma,I_{0},R,L$) correspond to the conductor phase (top) and the insulator phase (bottom). In top panels, for $Q_{+}^{2}+Q_{-}^{2}=1$ the core energy is not so large and CI is almost invisible. For $Q_{+}^{2}+Q_{-}^{2}=9$ we see the incoherence and the dissipation in the core region, signifying the CI. The conductor phase allows the proliferation of vortices but only those with $|Q_{\pm}|\leq 1$ are stable. In the bottom panels, both vortices have almost dissipated away due to EI, which starts from discrete poles near the boundary.555As a rule, it follows the sequence (51) found in the Appendix C from the pole structure of the propagator, though some of the steps can be absent, e.g. for a single Gaussian vortex there is no $\mathbb{C}_{2}$ stage. Indeed, the insulator phase has no free vortices, no matter what the charge. In Fig. 5 we see no instability even for a high-charge vortex in the perfect conductor phase (top), whereas the frustrated insulator phase (bottom) shows strong EI for the like-charged vortex $(3,3)$ since this fixed point has $g^{\prime}_{*}>0$, but the $(3,-3)$ vortex is stable. Notice that we could not expect CI for this case since the sum $Q_{+}^{2}+Q_{-}^{2}=9$ is the same in both cases – if for $Q_{-}=-Q_{+}$ the vortex has no CI, then for $Q_{-}=Q_{+}$ it cannot have it either (since the value $Q_{+}^{2}+Q_{-}^{2}$ is the same).

We have thus seen what patterns to expect from CI and EI, and also what kind of stable vortices to expect in different phases: the perfect conductor phase allows free proliferation of vortices of any charge, the conductor phase allows only single-quantum vortices (or vortices with sufficiently low $Q_{+}^{2}+Q_{-}^{2}$) while others dissipate from CI, the frustrated insulator supports the vortices with favorable charges (or favorable charge distribution in multiple-vortex systems) while others disintegrate from EI, and the insulator phase supports no vortices – they all dissipate from CI or EI, whichever settles first (depending on the vortex charges).

The case rich with analogies with condensed matter systems is the square vortex lattice on the background photonic square lattice, Fig. 6. Here we can also appreciate the transport processes. The photonic lattice is coincident with the beam lattice and equal in intensity, so $\Gamma(I_{0}+I_{x})=2\Gamma I_{0}$. In the perfect conductor phase (a) the vortices are stable and coherent and keep the uniform lattice structure. In the conductor phase (b) the CI is visible but the lattice structure survives. The bottom panels show the non-conducting phases, frustrated insulator (c) and insulator (d). Insulator looses both lattice periodicity and the Gaussian profile of the vortices but the frustrated insulator keeps regular structure: from EI the intensity is inverted and the resulting lattice is dual to the original one (compare (c) to (a)). The phase patterns $\theta_{F}(x,y;z=L/2)$ and $\theta_{F}(x,z;y=320\mu\mathrm{m})$ for the perfect conductor (top) and the frustrated insulator phase (bottom) are shown in the Fig. 7. Here we see the vortex charge transport mechanism in a PC: the vortices are connected in the sense that the phase $\theta_{F}$ is coherently traveling from one vortex to the next. In the FI phase, the phase is initially frozen along the $z$-axis, until the transport starts at some $z\approx L/2$.

It may be instructive to take a closer look at the lattice dynamics of the most interesting phase: the frustrated insulator. In Fig. 8 we inspect square lattices on the photonic lattice background for several charges of the form $(Q_{+}=3,Q_{-})$. The first row shows how the vortices loose stability and develop CI as the total square of the charge grows (from (a) to (c)). The panels ((d)-(i)) show how the $g^{\prime}$ coupling favors the opposite sign of $Q_{+}$ and $Q_{-}$ and how the optimal configuration is found for $Q_{-}=-3$. This is easily seen by minimizing the free energy over $Q_{-}$: it leads to the conclusion that the forward-backward coupling favors the ”antiferromagnetic” ordering in the sense that $Q_{+}+Q_{-}=0$.

Finally, it is interesting to see how the FI phase at high intensities and coupling strengths contains a seed of translation symmetry breaking which will become important in the presence of disorder. In Figs. 9-10 we give intensity and phase transverse profiles across the PC-FI transition and deep into the FI phase at large couplings. The intensity maps show the familiar inverse square lattice but the phase maps show stripe-like ordering, i.e. translation symmetry breaking along one direction in Figs. 10(c) and 10(d) – horizontal and vertical lines with a repeating constant value of the phase $\theta_{F}$ on all lattice cells along the line. This is a new instability, distinct from CI and EI. We cannot easily derive this instability from the perturbation theory in the Appendix C as it is a collective phenomenon and cannot be understood from a single beam.

IV The system with disorder

Consider now the same system in the presence of quenched disorder. This is a physically realistic situation: the disorder corresponds to the holes in the photonic lattice which are caused by the defects in the material. The defects are in fixed positions, i.e. they are quenched, whereas the beam is dynamical and can fluctuate. Now $I_{x}(\mathbf{r})\to I_{x}(\mathbf{r})+I_{h}(\mathbf{r})$, i. e., the quenched random part $I_{h}(\mathbf{r})$ is superimposed to the regular background (whose intensity is $I_{x}$). The disorder is given by some probability distribution, assuming no correlations between defects at different places. As in the disorder-free case, the lattice is static and ”hard”, i. e., does not backreact due to the presence of the beams. One should however bear in mind that the backreaction on the background lattice can sometimes be important as disregarding it violates the conservation of the angular momentum gaussrotpet . Disregarding the backreaction becomes exact when $I_{x}+I_{h}\gg|\Psi|^{2}$, i.e. when the background irradiation is much stronger than the propagating beams.

To treat the disorder we use the well-known replica formalism cardy . For vortex-free configurations, typical experimental values of the parameters suggest that the influence of disorder is small prlpet ; prapet ; vortlat2pet . However, the influence of disorder becomes dramatic when vortices are present. This is expected, since holes in the lattice can change the topology of the phase field $\theta_{\pm}$ (the phase now must wind around the holes). Our equations of motion are still given by the Lagrangian (4), but with $I_{x}\mapsto I_{x}+I_{h}$. In our analytical calculations, we assume that a defect in the photonic lattice changes the lattice intensity from $I_{x}$ to $I_{x}+I_{h}$, with Gaussian distribution of ”holes” in $I_{h}$, which translates to the approximately Gaussian distribution of the couplings $g,g^{\prime},g_{0},g_{1}$. In the numerics however, we do a further simplification and model the defects in a discrete way, i.e. at a given spot either there is a lattice cell of intensity $I_{1}$ (with probability $h$), or there isn’t (the intensity is zero, with probability $1-h$). This corresponds to $I_{x}=I_{1}/2,I_{h}=\pm I_{1}/2$ so the disorder is discrete. Due to the central limit theorem, we expect that the Gaussian analytics should be applicable to our numerics.

IV.1 The replica formalism at the mean-field level

To study the system with quenched disorder in the photonic lattice, we need to perform the replica calculation of the free energy of the vortex Hamiltonian (7). We refer the reader to the literature parisi ; spinglass for an in-depth explanation of the replica trick. In short, one needs to average over the various realizations of the disorder prior to calculating the partition function, i.e. prior to averaging over the dynamical degrees of freedom (vortices in our case). This means that we need to perform the disorder-average of the free energy, i.e. the logarithm of the original partition function $-\log\mathcal{Z}$, and not the partition function $\mathcal{Z}$ itself. The final twist is the identity $\log\mathcal{Z}=\lim_{n\to 0}\left(\mathcal{Z}^{n}-1\right)/n$: we study the Hamiltonian consisting of $n$ copies (replicas) of the original system and then carefully take the $n\to 0$ limit.666Care is needed as the $n\to 0$ limit does not in general commute with the thermodynamic limit. The partition function of the replicated Hamiltonian reads

[TABLE]

where $Q^{(\mu)}$ are the vortex charges in the $\mu$-th replica of the system. In the original Hamiltonian (7), the disorder turns the interaction constants into quenched random quantities $g_{ij},g^{\prime}_{ij},g_{0;ij},g_{1;ij}$, so we can compactly write our interaction term as

[TABLE]

with $J_{ij}^{++}=J_{ij}^{--}=g_{ij}(1-\delta_{ij})\log r_{ij}+g_{0}\delta_{ij}$, $J_{ij}^{+-}=J_{ij}^{-+}=g^{\prime}_{ij}(1-\delta_{ij})\log r_{ij}+g_{1}\delta_{ij}$. Now we again make the mean-field approximation for the long-ranged logarithmic interaction. Similar to the clean case, for $i\neq j$ we approximate $g\log r_{ij}\sim g^{\prime}\log r_{ij}\sim\log\Lambda$, knowing that $g,g^{\prime}\sim 1$ and assuming that average inter-vortex distance is of the same order of magnitude as the system size $\Lambda$, and for the core energy we likewise get $g_{0},g_{1}\sim\log a/\epsilon\sim-\log\epsilon\sim\log\Lambda$. The result is that all terms in $J_{ij}^{\alpha\beta}$, both for $i\neq j$ and $i=j$, are on average of the order $\log\Lambda\gg 1$, and the mean-field approach is justified. We will sometimes denote the $2\times 2$ matrices in the flavor space by hats (e.g. $\hat{J}=J^{\alpha\beta}$).

The final Hamiltonian (16) has the form of the random-coupling and random-field Ising-like model: random couplings stem from the stochasticity of $J_{ij}$ values and random field from the fact that $\langle J_{ij}\rangle\neq 0$ introduces terms linear in $Q_{i\alpha}$, i.e. an effective external field coupling to the ”spins”. We have arrived at this model through three steps of simplification: our microscopic model is a type of the XY-glass model (Cardy-Ostlund model cardyostlund ), a well-known toy model for disorder. At this stage our model is similar to the work of conti1 ; conti2 , only with two components instead of one. Then we have written the effective vortex Hamiltonian with Coulomb-like interaction, disregarding the topologically trivial configurations. This is a rather extreme approximation but a necessary one as it is very complicated to consider the full model with vortices. Finally, we have approximated the logarithmic potential with a constant all-to-all vortex coupling. Such an approximation (essentially the infinite dimension limit) is frequently taken and lies at the heart of the solvable Sherington-Kirkpatrick Ising random coupling model parisi . Our case differs from the Sherington-Kirkpatrick model as it (i) has also a random field (ii) has two flavors (iii) has the Ising spins taking arbitrary integer values. From the random XY-model it differs by (i) and (ii) above, and also by considering only vortices and no non-topological spin configurations. The additional phases we get in comparison to conti1 ; conti2 and its generalization in conti3 ; conti4 ; conti5 come from the interactions between the forward and backward flavor. But bearing in mind the drastic approximations we take we stress that we cannot aspire to solve neither the XY-model nor the resulting Ising-like model in any rigorous way (certainly not at the level of rigor of mathematical physics). We merely try to obtain a crude understanding.

The Gaussian distribution of defects reads $p(J^{\alpha\beta}_{ij})=\exp\left(-\left(J_{ij}^{\alpha\beta}-J^{\alpha\beta}_{0}\right)\left(\hat{\sigma}^{-2}\right)_{\alpha\beta}\left(J^{\alpha\beta}_{ij}-J^{\alpha\beta}_{0}\right)\right)$, where the second moments are contained in the matrix $\sigma_{\alpha\beta}$, with $\sigma_{+-}=\sigma_{-+}$. In this case we get the replicated partition function

[TABLE]

We can now integrate out the couplings $J^{\alpha\beta}_{ij}$ in (17) and get

[TABLE]

Integrating out the disorder has generated the non-local quartic term proportional to the elements of $\sigma_{\alpha\beta}^{2}$. The additional scale given by the average disorder concentration means we cannot scale out $\beta=L$ anymore, and it becomes an additional independent parameter. The partition function can be rewritten in the following way, usual in the spin-glass literature cardy ; spinglass . We can introduce the non-local order parameter fields

[TABLE]

which have the meaning of overlap between different metastable states. The rest is just algebra, although rather tedious: one rewrites the Hamiltonian in terms of new order parameters, and then one can solve the saddle point equations for $p_{\alpha}$ and $q_{\alpha\beta}$, or do an RG analysis. The calculation is found in Appendix G.

The mean-field analysis yields six phases:

One phase violates both the replica symmetry and the flavor symmetry, breaking it down to identity. We dub this phase vortex charge density wave (CDW), as it implies spatial modulation of the vortex charge, leading to nonzero net charge density $\sum_{i}Q_{i\alpha}^{(\mu)}$ in some parts of the system even if the boundary conditions are electrically neutral (the total net charge density must still be zero due to charge conservation). Vortices take their charges from $\mathbb{Z}\otimes\mathbb{Z}$. 2. 2.

The second phase violates the replica symmetry in both flavors and reduces the flavor symmetry but does not break it down to identity. Instead, it reduces it to the diagonal subgroup $U(1)_{F}\otimes U(1)_{B}\to U(1)_{d}$, so it has nonzero density of the vortex charge in a given replica $\sum_{i}Q_{i+}^{(\mu)}=-\sum_{i}Q_{i-}^{(\mu)}$. Again, the charge density is locally nonzero but now with an additional constraint resulting in frustration (multiple equivalent free energy minima!). This is thus the dirty equivalent of the frustrated insulator phase and we dub it vortex glass, as it has long-range correlations (because of the logarithmic interactions between charged areas), does not break spatial symmetry and exhibits frustration; its charges are from $\pi_{1}\left(U\left(1\right)_{d}\right)=\mathbb{Z}$. 3. 3.

The remaining phases have no nonzero vortex charge density fluctuation, and are similar to the phases in the clean system. Vortex perfect conductor violates the replica symmetry of all three fields $q^{++},q^{--},q^{+-}$ and allows free proliferation of vortices with charges $(Q_{+},Q_{-})\in\mathbb{Z}\otimes\mathbb{Z}$. 4. 4.

Frustrated vortex insulator preserves the replica symmetry of $q^{\pm\pm}$ but has non-zero value, with broken replica symmetry, of the mixed $q^{+-}$ field, which gives $U(1)_{d}$ vortices, with charges $Q_{+}=-Q_{-}\in\mathbb{Z}$. 5. 5.

Vortex conductor preserves the replica symmetry of the mixed $q^{+-}$ order parameter but violates it in $q^{\pm\pm}$, resulting in the proliferation of single-flavor vortices with $\mathbb{Z}$ charge. 6. 6.

Vortex insulator fully preserves the replica symmetry, all order parameters are zero and vortices cannot proliferate. RG analysis will show that insulator surivives only at zero disorder, otherwise it generically becomes CDW.

The phase diagram (given in Fig. 11 in the next subsection) now contains six phases (only five are visible for the parameters chosen in the figure): CDW, insulator, FI, conductor, PC and the glassy phase. The insulator phase is now of measure zero in the $(g,g^{\prime},\sigma^{2})$ plane, existing only for the points at $\sigma^{2}=0$; for generic nonzero values we have a CDW. For simplicity, we have plotted the phase diagram for $\sigma_{++}^{2}=\sigma_{--}^{2}=\sigma_{+-}^{2}\equiv\sigma^{2}$.

IV.2 RG analysis and the phase diagram

To study the RG flow, we can start from the replicated partition function (18), inserting the definition of the couplings $J_{ij}^{\alpha\beta}$ and keeping the vortex charges $Q^{(\mu)}_{i\alpha}$ as the degrees of freedom (without introducing the quantities $p^{\mu}_{\alpha},q^{(\mu\nu)}_{\alpha\beta}$). The basic idea is the same: we consider the fluctuation $\delta\left(\bar{\mathcal{Z}}^{n}\right)$ upon the creation of a vortex pair at $\mathbf{r}_{1,2}$ with charges $\vec{q}^{(\mu)}_{1},-\vec{q}^{(\mu)}_{2}$, in the background of the vortices $\vec{Q}^{(\nu)}_{1,2}$ at positions $\mathbf{R}_{1,2}$. Likewise, we introduce the fugacity parameter $y^{(\mu)}$ to account for the vortex core energy. However, this problem is much harder than the clean problem and one has to resort to many approximations to perform the calculation. In its most general form, the problem is still open, in the sense that all known solutions suppose a certain form of replica symmetry breaking or truncate the RG equations spinglass . The RG analysis is thus less useful in the disordered case but at least the numerical integration of the flow equations is supposed to give a more precise rendering of the phase diagram compared to the mean field theory. We again describe the calculation in Appendix G and jump to the results.

The fixed point of the flow equations lies either at infinite $y$ or at $y=0$ like in the clean case. This is again controlled by the the equation for $\partial y/\partial\ell$ but now depending on the combination $g+g^{\prime}+\beta^{2}\sigma^{2}$ instead of $g+g^{\prime}$ in the clean case (for simplicity, we consider the case where $\sigma_{\alpha\beta}^{2}$ are all equal). The following cases appear:

When the fugacity flows toward infinity, we reproduce the phases and the fixed point values $(g,g^{\prime},\sigma^{2})$ from the clean case: the PC flows toward $(0,0,0)$, the FI toward $(g_{*},g^{\prime}_{*},0)$ and the conductor toward $(g_{*},g^{\prime}_{**},0)$ with $g_{*}\to-\infty,g^{\prime}_{*}\to-\infty,g^{\prime}_{**}\to\infty$. Notice that all these phases flow to $\sigma^{2}=0$, i.e. disorder is irrelevant. 2. 2.

When the fixed point lies at $y=0$, one possibility is that all parameters ($g,g^{\prime},\sigma^{2}$) flow toward some non-universal nonzero values. The attraction region of this point is the CDW phase: the disorder term stays finite as well as the couplings. In particular, the points on the half-plane $g+g^{\prime}>0,\sigma^{2}=0$ stay at $\sigma^{2}=0$ (with constant coupling values) and this is the insulator phase from the clean case. Notice that $\sigma^{2}>0$ now, i.e. disorder is relevant. For $\sigma^{2}<1$, this are the only fixed points when $y=0$. 3. 3.

However, for sufficiently strong disorder ($\sigma^{2}>1$), there is a new line of fixed points at $y=0$ with a finite attraction region, corresponding to a new phase. For $\beta>1$, the right-hand side of the second RG equation in (93) has a zero at nonzero $g^{\prime}$ and there are trajectories flowing toward $(y,g,g^{\prime},\sigma^{2})=\left(0,g,g^{\prime}(g),\sigma^{2}(g)\right)$ and not toward an arbitrary nonuniversal value of $\sigma^{2}$. This is precisely the glass phase, where disorder is again relevant. At the lowest order, the relation between $g,g^{\prime},\sigma^{2}$ at the fixed point line is given by the relation $g+g^{\prime}+\beta^{2}\sigma^{2}=1$.

Now we have made contact between the mean-field classification of phases and the fixed points and regions of the RG flow. The flows in the $(g,g^{\prime})$ plane are given in the Fig. 11. The parameter space is four-dimensional so the phase structure is different at different disorder concentrations $\sigma^{2}$. In the (A) panel for $\sigma^{2}=0.4$, the phase structure is similar to the clean case – we see the same four phases except that insulator (no stable vortices) is replaced by the CDW phase with localized vortices. In the (B) panel, for $\sigma^{2}=1.2$, the CDW phase is replaced by another disordered phase, the glass-like regime. Importantly, the glass phase does not cross the $g^{\prime}=0$ axis, meaning that a single-flavor system even with disorder could not support a glass. We thus conjecture that the transition at $\sigma^{2}=1$ is of first order, as the change is the structure of the $(g,g^{\prime})$ phase diagram is discontinuous, and we do not see how this could happen if the first derivative $\partial\mathcal{F}/\partial\rho_{\pm}$ (the derivative of the free energy with respect to vortex charge density) is continuous. However, we have not checked the order of this transition by explicit calculation. The phase structure is further seen in the $\sigma^{2}-g^{\prime}$ diagram, where we see the glass phase emerge at some value of the disorder. This is discussed further in the next section, where we study the equivalent antiferromagnetic system (with the same structure of the phase diagram, Fig. 16).

IV.3 Geometry of patterns

The two previously considered mechanisms of instability – central instability and edge instability – remain active also in the presence of disorder. However, in the presence of disorder there is a third, inherently collective effect that we dub domain instability (DI). It follows from the fact that the self-focusing term $\Gamma E$ grows with intensity $I$: more illuminated regions react faster (Eqs. 1-2). In the presence of background lattice there will be regions of initially zero beam intensity $I_{0}$ where the regular lattice cells have some nonzero intensity $I_{x}$. Approximating $I=I_{0}+I_{x}\approx I_{x}=\mathrm{const.}$, our equations in the vicinity of the defect (hole) in the background lattice becomes the Schrödinger equation in a step potential (equal to $I_{x}$ in the regular parts of the photonic lattice, and equal to zero where a hole is found), so the $z$-dependent part of the solution is of the form $\sum_{k}e^{\imath\lambda_{k}z}$ and the eigenenergies along $z$ are gapped by the inverse length: $\lambda_{k}>1/L$. For small eigenenergies, the transmission coefficient is very low whereas for large energies it approaches unity. Thus for $1/L$ large (i.e., there are few $\lambda_{k}$’s which are larger than $1/L$) most of the intensity remains confined by the borders of the defect and the intensity does not spill but for small $1/L$ the beam profile is deformed by the ”spilling” into the hole regions. For vortices, there is an additional Coulomb interaction in the $x-y$ plane, meaning the effective potential is not piecewise constant anymore (even in the simplest approximation) but the qualitative conclusion remains: large $L$ brings global reshaping of the intensity profile.

The other phases are analogous to the ones in the clean case, though with a general trend that the presence of disorder decreases the stability of vortex patterns. The PC and FI phases are shown in Figs. 12-13. In this section we only look at the lattices, as the notion of disorder is inapplicable for a single beam. Consider first the patterns in the PC phase (Fig. 12). Compared to the clean case (Fig. 6(a)), the symmetry is much reduced, from $O(2)$ to $C_{4}$ but the vortices are conserved and the original lattice structure (outside the holes) is clearly visible. The FI (Fig. 13) shows mainly EI (and to a smaller extent CI), which together lead to the lattice inversion. The rule of thumb for differentiating the conductor and PC on one side from the CDW and FI on the other side is precisely the presence of the lattice inversion. The absence of the charge transport is best appreciated in the phase images: the charge pins to the defects and localizes toward the end of the crystal (i. e., for $z$ near $L$). Only near the edges we see high vorticity, somewhat analogous to topological insulators, which only have nonzero conductivity along the edges of the system.

The CDW versus the glass phase is given in Fig. 14. The charge density wave (a, c, $L=240\mu\mathrm{m}$) exhibits the diffusion of intensity due to DI, and the vortex beams are in general asymmetric and not clearly delineated. In (b, d), where $L=120\mu\mathrm{m}$ with all other parameters the same, there is a clear border between defects and the regular parts of the lattice and the intensity is concentrated in the vortex cores. We give also the vortex charge density map in (c, d) in addition to the intensity maps in (a, b) as the charge density shows why the CDW is insulating: even though individual beams diffuse and smear out in intensity, the regions of nonzero vortex charge are disjoint and no global conduction can occur. Glass is divided into ordered domains in intensity but the vortex charges form a connected network which supports transport. This is analogous to the percolation transition in a disordered Ising model parisi1 ; parisi2 and we may expect that the CDW-glass transition follows the same scaling laws near the critical point. However, we have not checked this explicitly and we leave it for further work.

V The condensed matter analogy: collinear doped Heisenberg antiferromagnet

The two-beam photorefractive system can serve as a good model for quantum magnetic systems. The most obvious connection is to multi-component XY antiferromagnets (i.e., two-dimensional Heisenberg model): planar spins are nothing but complex scalars, and the vortex Hamiltonian remains identical ($\pi_{1}(SO(2))=\pi_{1}(U(1))=\mathbb{Z}$). The nonlinearity in the spin system is different and usually much simpler, but that typically does not influence the phase diagram (the symmetry structure remains the same). Such connection is so obvious it does not require further explanations. Our point is that the CP beams in a PR crystal can also describe more general magnetic systems in the presence of topological solutions described by homotopy groups different from $\mathbb{Z}$. In particular, we want to point out to a connection with a two-sublattice antiferromagnetic system which has some time ago enjoyed considerable popularity as a possible description of magnetic ordering in numerous planar strongly-coupled electron systems, including cuprate high-$T_{c}$ superconductors ft2 ; sachbook ; modelza . This is the collinear doped antiferromagnet defined on two sublattices. When coupled to a charge density wave (speaking about the usual, $U(1)$ electromagnetic charge) and a superconducting order parameter, it becomes a toy model of cuprate materials (one variant is given in modelza ). In the light of what we know today, the ability of this model to realistically describe the cuprate physics is quite questionable; but even so it is an interesting magnetic system on its own, and it was already found in jurms0 ; jurms1 to exhibit a spin-glass phase, though in a slightly different variant (in particular, with spiral instead of collinear ordering).

Let us formulate the model. While the material is a lattice on the microscopic level, here we are talking about an effective field theory model. The order parameter is the staggered magnetization

[TABLE]

where $\alpha\in\{1,2\}$ is the sublattice ”flavor” index (analogous to the $\alpha$-index for the $F$ and $B$ beam in the previous sections)777Sometimes we will denote the sublattices by $\pm$ instead of $1,2$ for compactness of notation. and each component $\mathbf{M}_{\alpha}$ is a three-component spin, describing the internal, i.e. spin degree of freedom (we label the spin axes as X,Y,Z). The total spin is thus the sum of the spins of the two components, and $\mathbf{n}$ is the modulation vector. The modulation gives rows of alternating staggered magnetization in opposite directions as in Fig. 15(a). This stands in contrast with the spiral order, where the modulation vectors become $\mathbf{n}_{\alpha}$, i.e. differ for the two sublattices, and are themselves space-dependent jurms0 . The ordered phase of the collinear system has the nonzero expectation value of the staggered magnetization along one direction, which can be chosen as the Z-axis (”easy axis”), where the spin fluctuations about the easy axis remain massless, and the symmetry is broken from $O(3)$ to $O(3)/O(2)$. The spiral order, on the other hand, breaks the symmetry down to identity, as the order parameter is a dreibein jurms0 .

The symmetry conditions (isotropy in absence of external magnetic field) determine the Hamiltonian up to fourth order, as discussed in modelza :

[TABLE]

The antiferromagnetic coupling is $g_{M}$, the spin stiffness is $c_{M}$ and the effective mass of spin wave excitations is $r$. The fourth-order coupling $u_{0}$ comes from the ”soft” implementation of the constraint $|\mathbf{M}_{\alpha}|=1$888One could also enforce the constraint exactly, through the nonlinear sigma model, as was done in jurms0 . While the leading term of the ”vortex” Hamiltonian would remain the same in that case, the amplitude fluctuations have different dynamics which influences some terms of the Hamiltonian and thus its RG flow (though probably not the very existence of the glass phase). and $v_{0}$ is the anisotropy between the two sublattices, justified by the microscopic physics ft2 ; modelza . The Hamiltonian can be transformed by rescaling $\tau$ and $x,y$, together with the couplings $u_{0}\mapsto u$ and $v\mapsto v_{0}$ to set $g_{M}=c_{M}=1$ so that the kinetic term becomes isotropic, giving

[TABLE]

where we have also rewritten the quartic terms for convenience. Without anisotropy, the energy of the system is a function of $|\mathbf{M}_{1}|^{2}+|\mathbf{M}_{2}|^{2}$ only and the symmetry group is the full $O(6)$. With $v\neq 0$, the symmetry is reduced to $O(3)_{1}\otimes O(3)_{2}$: the internal spin symmetry in each sublattice remains unbroken but the spatial rotation symmetry between the layers is broken down to just the discrete flip. Compare this to the $U(1)\otimes U(1)$ symmetry in the PR system: there, it is the internal phase symmetry that remains unbroken.

V.1 $\mathbb{Z}_{2}$ ”vortices”

Remembering that topological solitons are classified by homotopy groups, and that we work in a two-dimensional plane, the relevant group is again the first homotopy group, $\pi_{1}(O(3))=\mathbb{Z}_{2}$. For simplicity, we will call these excitations ”vortices”, bearing in mind that the only possible charges are $Q_{\alpha}=\pm 1$ and not all integers. A realization of the ”vortex” with $Q=1$ is shown in Fig. 15(b). Since the spins are three-dimensional (the figure shows the projection in the XY plane), it becomes clear that ”vortex” charge is only defined modulo $2$, i.e. it makes no sense to talk about charges $|Q|>1$. For example, winding around twice in the XY plane can be done along a closed line in the XYZ space which can be contracted to a point. That could not happen for the two-dimensional phase $U(1)$ precisely because there is no extra dimension. In Fig. 15(b) the ”vortex” is superimposed onto the regular configuration: it is recognizable as a contact point between two lines of alternating staggered magnetization. In. Fig. 15(c) we have subtracted the regular part and only the ”vortexing” spin pattern is shown: here we see the ”vortex” interpolates between two opposite spin orientations in two opposite directions in the plane.

Now let us derive the effective Hamiltonian of the ”vortices”. For the $\mathbb{Z}_{2}$ ”vortex”, a loop in real space is mapped onto a $\pi$-arc in the internal space, so the ”vortex” can be represented as

[TABLE]

giving (the matrices $\ell_{1,2,3}$ represent the $so(3)$ algebra):

[TABLE]

where $\mathbf{m}_{\alpha}$ is the magnetization amplitude, analogous to the beam amplitude $\psi_{\alpha}$ in the optical system. The leading-order, non-interacting term in (22) gives for the energy of a single ”vortex” of charge $\vec{Q}$:

[TABLE]

which is in fact independent of the sign of $\vec{Q}$ (as could be expected, as it is in general proportional to $\vec{Q}\cdot\vec{Q}$ which is a constant for parity ”vortices”). The ”vortex” singles out an easy axis (Z-axis) around which the staggered magnetization winds ($\phi$ being the winding angle). This allows one to introduce $\mathbf{m}_{\alpha\perp}\equiv\left(m_{X\alpha},m_{Y\alpha},0\right)$. A ”vortex” pair with charges $\vec{Q}_{i}$ and $\vec{Q}_{j}$ has the binding energy

[TABLE]

Now we should integrate out the amplitude fluctuations as we did in Appendix D for the CP beams. This again leads to the coupling between different flavors, giving a ”vortex” Hamiltonian analogous to (7):

[TABLE]

Two obvious differences with respect to the optical system are (i) the charges are now limited to the values $\pm 1$ (ii) there is a term linear in charge density, which acts as a chemical potential. The latter arises from the coupling of the three-dimensional spin waves (i.e., the topologically trivial excitations of the amplitude $\mathbf{m}_{\alpha}$) to the ”vortices”. Remember that in the CP system, the amplitude fluctuations also couple to the vortices, but there is no third, $\mathrm{Z}$-axis of the order parameter so no linear term appears. The microscopic expressions for the effective parameters $g,g^{\prime},\mu_{\alpha}$ read:

[TABLE]

assuming $m_{1\perp}=m_{2\perp}\equiv m_{\perp}$. Now the RG calculation is similar to the optical case but the nonzero chemical potential introduces two differences. First, there is obviously the additional term proportional to the total charge of the virtual pair of ”vortices”, $\mu_{\alpha}(q_{1\alpha}+q_{2\alpha})$. Second, there is no charge conservation as the expectation value of the total ”vortex” charge is now $\langle\vec{Q}\rangle=\partial\mathcal{F}/\partial\vec{\mu}\neq 0$. Thus we need to take into account not only the fluctuations with zero net charge (virtual ”vortex” pairs with charges $\vec{q}_{1}\equiv\vec{q}$ and $\vec{q}_{2}\equiv-\vec{q}$) but also the situations with arbitrary pairs $\vec{q}_{1},\vec{q}_{2}$.999In the CP beam system, the total vortex charge can be nonzero if the boundary conditions at $z=0,L$ have nonzero total vorticity. But there we had no bulk chemical potential so the total vorticity in the crystal could not change during the propagation along $z$. Here, we have a bulk term in the Hamiltonian which violates charge conservation. This modifies the variation of the partition function from (12-13) to:

[TABLE]

where we have introduced $2\vec{q}\equiv\vec{q}_{1}-\vec{q}_{2},\vec{q}_{0}\equiv\vec{q}_{1}+\vec{q}_{2}$ and $\delta\mathbf{R}_{1,2}\equiv\mathbf{R}_{1,2}-\mathbf{r}$. The mixed term which includes both $\vec{q}$ and $\vec{q}_{0}$ vanishes due to isotropy. Matching the terms in the resulting expression with the original Hamiltonian, we find the recursion relations:

[TABLE]

Crucially, the chemical potential does not run which could be guessed from dimensional analysis (it couples to dimensionless charge). This is the same system as (14) up to the trivial rescaling of the coupling constants and the shift of the critical line $g+g^{\prime}=1$ in the PR system to the line $g+g^{\prime}+\mu_{+}+\mu_{-}=1$. It becomes obvious that the phase diagrams are equivalent and can be mapped onto each other.

V.2 Influence of disorder

The disorder in a doped antiferromagnet comes from electrically neutral metallic grains quenched in the bipartite lattice. Being metallic and neutral, they are naturally modeled as magnetic dipoles $\mathbf{X}$ quenched in the bipartite lattice. This picture stems from the microscopic considerations in jurms3 . We again assume the Gaussian distribution of the disorder as $p(X)\propto\exp(-|\mathbf{X}|^{2}/2\sigma_{X}^{2})$. The disorder dipoles are one and the same for both sublattices, so $\mathbf{X}$ has no flavor (sublattice) index. The minimal coupling of the dipoles to the lattice spins $\partial_{i}\mapsto\partial_{i}-\imath\hat{\ell}_{i}X_{i}$ gives

[TABLE]

Now the replica calculation requires the multiplication of the $M$ field into $n$ copies and performing the Gaussian integral over the disorder. The initial distribution of the disorder $p(X)$ gives rise to two independent Gaussian distributions: for the couplings $J_{ij}^{\alpha\beta}$ with dispersion matrix $\sigma^{2}_{\alpha\beta}$, and for the chemical potential $\mu_{i}^{\alpha}$ with the dispersion vector $\xi^{2}_{\alpha}$. The resulting Hamiltonian is

[TABLE]

where we have disregarded the subleading logarithmic term ($\sim\log|\mathbf{M}_{\alpha}^{(\mu)}|$). Now making use of the representation (23) and plugging it in into (33) gives the disordered ”vortex” Hamiltonian

[TABLE]

Of course, we could have arrived at the same effective action starting from the ”vortex” Hamiltonian (27), taking the infinite range approximation and identifying $J_{ij}^{\alpha\alpha}=g_{ij}\log r_{ij}$ and similarly for other components of $J_{ij}^{\alpha\beta}$ as we demonstrated for the PR system. The final result has to be same at leading order.

The next step is to rewrite the Hamiltonian in terms of the order parameters $p_{\alpha}^{(\mu)},q_{\alpha\beta}^{(\mu\nu)}$ defined in (19). Compared to the effective action for the photonic lattice with disorder in Eq. (78), there are two extra terms in the resulting action $S_{\mathrm{eff}}$: one is proportional to the dispersion $\xi^{2}$ and the other to the mean chemical potential $\vec{\mu}_{0}$. The former term just introduces the shift $J_{0}^{\alpha\beta}\mapsto J_{0}^{\alpha\beta}-\sigma^{2}/2\beta$ and the latter term, linear in the ”vortex” charges and proportional to the chemical potential, introduces solutions with nonzero net ”vortex” charge density. Looking back at the results of the saddle-point calculation in Eqs. (19,88), this tells us that the relation between the phase diagrams is the following. The phases with no net ”vortex” charge density – insulator, conductor, frustrated insulator and perfect conductor – remain the same as in the PR system, since both the average coupling value $J_{0}^{\alpha\beta}$ (which gets shifted) and the term proportional to the chemical potential $\mu_{\alpha}$ couple only to $\vec{p}^{(\mu)}$. For brevity, denote $J_{0}^{\pm\pm}\equiv J_{0}^{\pm}$ and notice that $J_{0}^{-+}=J_{0}^{+-}$. The structure of phases with nonzero $\vec{p}^{(\mu)}$ depends on the zeros of the saddle point equation

[TABLE]

analogous to (87), where the 1-step replica symmetry breaking implies $p_{(\mu)}^{\pm}=\left(p^{\pm},\ldots,p^{\pm}\right)$. Now the equation is cubic and the structure of solutions is different from (88). We could not find the solution in the closed form but it is clear that a pair of cubic equations will have either a single solution $(p^{+},p^{-})$ or nine combinations $(p^{+},p^{-})$, not necessarily all different. Numerical analysis of (35) reveals only two inequivalent solutions, analogous to (88), i.e. one of them has a single free energy minimum, the other one a pair of degenerate minima. Therefore, we again have two disordered solutions, one of which is glassy (frustrated).

Now we can write down also the RG equations for the effective action (34). In this calculation, we put $\xi_{\alpha}^{2}=\sigma_{\alpha\beta}^{2}\equiv\sigma^{2}$ for simplicity. Following the same logic as earlier, the equations are found to be:101010For the most general case of different and non-scalar $\sigma^{2}_{\alpha\beta}$ and $\xi^{2}_{\alpha}$ the flow equations for them complicate significantly and we will not consider them.

[TABLE]

Like in the clean case, the chemical potential is irrelevant and the solutions for fixed point are the same as for the PR beams, including the spin-glass fixed point. We conclude that the phase structure of the optical system is repeated in strongly correlated doped antiferromagnets, which also exhibit the spin-glass phase and have the phase diagram sketched in Fig. 16. In this context it is more interesting to plot the phase diagram in the $\sigma^{2}-1/g^{\prime}$ plane, mimicking the $x-T$ phase diagram of quantum critical systems sachbook (remember that the coupling constants $g,g^{\prime}$ behave roughly as inverse temperature in XY-like models). Bear in mind that all phases shown are about ”vortex” dynamics, i.e. one should not compare Fig. 16 to the textbook phase diagram of high-temperature superconductors, which accounts also for the charge or stripe order and the superconducting order. All ”vortex” phases would be located inside the pseudogap regime of the superconductor, where various exotic orders can coexist (assuming, of course, that our model is an adequate approximation of the magnetic order in a cuprate or similar material, which is a complex question). Crucially, the spin-glass phase (blue curves) flows toward finite disorder $\sigma^{2}$, whereas the remaining two phases end up at zero disorder, either at infinite $1/g^{\prime}$ (PC, red flows), or at zero $1/g^{\prime}$ (conductor, green flows). The RG flows in the conductor phase are almost invisible in the figure, as the flows are much slower than in the remaining two phases.

V.2.1 Discussion

Early papers which found and explored the spin-glass phase in a very similar model are jurms0 ; jurms1 ; jurms2 ; jurms4 . The main difference is that the papers cited consider the spiral (non-collinear) spin order. These works are all inspired by the cuprate materials, the most celebrated brand of high-temperature superconductors. While jurms1 ; jurms2 explore in detail the transport properties, we have no pretension neither to provide a realistic model of cuprates nor to explore in detail all the properties of the spin-glass phase. We are content to see that the PR system of $\mathbb{Z}$ vortices reproduces the phase structure of a certain kind of dirty Heisenberg antiferromagnets (with $O(3)$ spins and $\mathbb{Z}_{2}$ ”vortices”), besides the more obvious connection to systems which directly reproduce the $\mathbb{Z}$ vortices in multicomponent $U(1)$ systems like multicomponent Bose-Einstein condensates and type-$1.5$ superconductors.

VI Conclusions

We have investigated the light intensity patterns in a nonlinear optical system, a pair of counterpropagating laser beams in a photorefractive crystal. We have studied this system as a strongly interacting field theory and have focused mostly on the formation and dynamics of vortices. The vortices show a remarkable collective behavior and their patterns are naturally classified in the framework of statistical field theory: the effective action shows several different phases with appropriate order parameters, and the system is essentially an XY model with two flavors, i.e. two kinds of vortex charge, for the two beams. The interaction between the flavors is the central reason that the total energy of the Coulombic interactions between the vortices in general cannot be locally minimized at every point. In the presence of disorder, a phase with multiple free energy minima arises, where the absence of long-range order is complemented by the local islands of ordered vortex structure, and which resembles spin glasses.

The phase diagram is simple in terms of the effective parameters – vortex coupling constants – and quite complex when expressed in terms of the experimentally controllable quantities – the intensity of the laser beams, the intensity of the background photonic lattice and the properties of the photorefractive crystal (the last is not controllable but can be estimated reasonably well korn ). The lesson is that the approach we adopt can save us from demanding numerical work if the space of original parameters is ”blindly” explored. Our phase diagrams can serve as a starting point for ”guided” numerical simulations, suggesting what phenomena one should specifically look for. So far the field-theoretical and statistical approach was not much used in nonlinear optics (important exceptions are primeri1 ; primeri2 ; fisher ; fisher2 ; fisher3 ; fisher4 ; conti1 ; conti2 ; conti3 ; conti4 ). We hope to stimulate work in this direction, which is promising also because of the potential of the photorefractive systems to serve as models of strongly correlated condensed matter systems. They make an excellent testing ground for various models because of the availability and relatively low cost of experiments.

In this work we have focused on the relation of the photorefractive counterpropagating system to the model of an $O(3)$ doped antiferromagnet with two sublattices. The authors of previous works on this model modelza ; jurms1 ; jurms2 ; jurms3 were motivated mainly by the ubiquitous problem of understanding the pseudogap phase in cuprate superconductors. The applicability of the model to this particular problem is still an open question – it may well be that cuprate physics goes far beyond. Nevertheless, it is an important quantum magnetic system in its own right and serves as an illustration on how one can simulate condensed matter systems in photorefractive optics.

Another field where vortices are found as solutions of a non-linear Schrödinger equation are cold atom systems and Bose-Einstein condensates bec . Notice however that Bose-Einstein condensates in optical traps are usually (but not always, see becmulti ) three-dimensional systems with vortex lines (rather than XY-type systems with point vortices) and our formalism would be more complicated there: in three spatial dimensions, vortices give rise to emergent gauge fields. The multi-component systems of this kind give rise to so-called type $1.5$ superconductors multi15 which are a natural goal of further study.

A more complete characterization of the glass-like phase is also left for further work. The reader will notice we have devoted very little attention to the correlation functions in various parameter regimes or the scaling properties of susceptibility, which should further corroborate the glassy character of the system. This is quite difficult in general but very exciting as it offers an opportunity to tune the parameters (e.g. disorder strength) freely in the optical system and study the glass-like phase and its dynamics.

Acknowledgments

We are grateful to Mariya Medvedyeva and Vladimir Juričić for careful reading of the manuscript. Work at the Institute of Physics is funded by Ministry of Education, Science and Technological Development, under grants no. OI171033 and OI171017. M. P. is also supported by the NPRP 7-665-1-125 project of the Qatar National Research Fund (a member of the Qatar Foundation).

Appendix A Numerical algorithm

In order to solve numerically the system (1-2) we employ a variation on the method of belicalg0 ; belicalg1 . The method does not make use of any analytical ansatz: it is an ab initio numerical procedure which integrates the equations of motion. The system has four independent variables: the transverse coordinates $(x,y)$, the longitudinal coordinate (”formal time”) $z$ and the (physical) time $t$. That means we have essentially three nested loops: (i) at every $z$-slice we integrate the transverse Laplacian and the interaction terms (ii) when done with this for the whole $z$-axis we advance the time $t$ (iii) we repeat the whole procedure until reaching some time $t_{f}$ which certainly should be much longer than the relaxation time $\tau$.

The important point is the very different nature of the initial and boundary conditions for various coordinates. The boundary conditions in the $(x,y)$ plane, i.e. at the crystal edge are not crucial: we have either just one or a few Gaussian beams whose intensity drops exponentially away from the center and is practically zero at the crystal edge, or we have a large lattice consisting of many (of the order of $50-100$) Gaussian beams so the edge effects only affect a small portion of the whole lattice. Therefore, imposing periodic boundary conditions (stemming naturally from the integration in Fourier space, see the next paragraph) are perfectly satisfying. Crucially, however, the CP geometry means that $F(t;z=0;x,y)=F_{0}(x,y)$ and $B(t;z=L;x,y)=B_{0}(x,y)$ are given functions, fixed for all times. We thus have a two-point boundary value problem along $z$ and have to iterate the $z$-integration several times until we reach the right solution. Finally, the initial condition for the relaxation equation (2) is that the crystal is initially at equilibrium, meaning that $E(t=0)=-I_{x}/(1+I_{x})$; specifically, for zero background lattice $E(t=0)=0$.

The algorithm now has the following structure:

The innermost loop integrates in the $x-y$ plane. This is a Poison-type (elliptic) equation, thus we employ the operator-split method, integrating the Laplacian operator in the Fourier space and the interaction term (the $EF$ and $EB$ terms) in real space, in the second-order leapfrog scheme. Thus, at every time instant $t_{i}=i\Delta t$, we start from $z=0$ where we set the condition $F(i\Delta t;z=0;x,y)=F_{0}(x,y)$, divide the $z$-axis into $N$ steps of size $\Delta z=L/N$, and at every slice $z=j\Delta z$ we perform the frog’s leap: we do the Fast Fourier Transform (FFT) to turn the $(x,y)$-dependence into $(q_{x},q_{y})$-dependence,111111We denote the fields in Fourier space with a tilde, e.g. $\tilde{F}$. then we advance the Laplacian for $\Delta z/2$ as $F(i\Delta t;j\Delta z;\mathbf{q})\equiv\tilde{F}_{i,j}^{(0)}\mapsto\tilde{F}_{i,j}^{(1)}=\exp\left(-\imath q^{2}\Delta z/2\right)\tilde{F}_{i,j}^{(0)}$, then we do the inverse FFT and advance the interaction in real space as $F_{i,j}^{(2)}=\exp\left(\imath\Gamma E\left(i\Delta t;j\Delta z;x,y\right)\right)F_{i,j}^{(1)}$. Finally we do the FFT again and advance the Laplacian for the remaining half-step, $\tilde{F}_{i,j+1}=\exp\left(-\imath q^{2}\Delta z/2\right)\tilde{F}_{i,j}^{(2)}$. Once we reach $j=N$ the integration goes backward, along the same lines, updating now the $B$-field (starting from $B_{0}(x,y)$), where all signs in the exponents of the above formulas are to be reversed. When we reach $z=0$ again, we are done. In this loop we use the field $E_{1,j}$ as already known for all $j$. 2. 2.

The above loop will in general produce results inconsistent with the charge field $E_{i,j}$ because the equation for $E$ couples $F$ and $B$ and we have ignored that by integrating the two fields one after the other instead of simultaneously. This is, of course, commonplace in two-point boundary value problems: either only one boundary condition can be imposed exactly and the other is shot for or, as in our case, both are imposed exactly but at the cost of the solution being inconsistent with the equations, so we have to iterate the system to arrive at the correct solution everywhere. The second loop thus iterates the first loop $A$ times, at each step updating the charge field as $E_{i,j}^{(a-1)}\mapsto E_{i,j}^{(a)}=E_{i-1,j}-\tau\left(E_{i,j}^{(a-1)}+I_{i,j}^{(a-1)}/\left(1+I_{i,j}^{(a-1)}\right)\right)/(1+I_{i,j}^{(a-1)})$. The number of iterations $A$ is not fixed: we stop iterations when the intensity pattern stabilizes, $\sum_{j}\sum_{x,y}\left(I_{i,j}^{(a)}-I_{i,j}^{(a-1)}\right)<\epsilon$, for some tolerance $\epsilon$. Here, $I_{i,j}$ refers to total intensity, i.e. $|F|^{2}+|B|^{2}+I_{x}$. 3. 3.

Finally, the outermost loop integrates in time $t$, from $t=0$, with the initial condition $E(t=0)=-I_{x}/(1+I_{x})$ given above. The integration time $t_{f}$ is divided into $M=t_{f}/\Delta t$ intervals, and at the end of each step we update $\left(F_{i,j},B_{i,j},E_{i,j}\right)\mapsto\left(F_{i+1,j},B_{i+1,j},E_{i+1,j}\right)$. Only the charge field is directly integrated (as written above), in the first-order, Euler scheme. The beam envelopes depend on time only parametrically, through $E(t)$, and they evolve by using an updated $E_{i,j}$ in the first two loops at every time step.

This procedure is very close to that in belicalg0 ; the main difference is that we use a second-order (leapfrog) scheme while on the other hand our time integration is of the lowest, linear order instead of second order as in belicalg0 .

Appendix B Time-dependent perturbation theory and the existence of equilibrium configurations

B.1 Stability analysis: fixed points and limit cycles

In this Appendix we consider the time evolution of the CP beams and show the existence of a stable equilibrium point with nonzero intensity. This means that the system reaches stationary state for long times, justifying the basic assumption of the paper that one can study the vortex configurations within equilibrium statistical mechanics. Not all patterns are stable: depending on the boundary conditions and parameter values, the system may or may not have a stable equilibrium, and non-equilibrium solutions in photorefractive optics are well-known fisher ; gaussrotpet . For our purposes, however, it is enough to identify the region of parameter space where the equilibrium exists; other cases are not the topic of this paper.

The time evolution of the beams $\Psi_{\alpha}$ and the charge field $E$ in $(k,q)$-space is obtained by differentiating the equations (1) with respect to time and plugging in $\partial E/\partial t$ from the relaxation equation (2):

[TABLE]

This system has three equilibrium points. One is the ”[math]” point,

[TABLE]

and the remaining two are related by a discrete symmetry $\Psi_{\pm}\mapsto\Psi_{\mp}$, so we denote them as ”$\pm$” points, the ”$+$” point being

[TABLE]

and the ”$-$” point has instead $\Psi_{+}=0$ and $\Psi_{-}=\sqrt{(E(1+I_{x})+I_{x})/(1+E)}\exp(\imath\phi_{-})$. Notice that the phase $\phi_{\pm}$ remains free to vary so this solution supports vortices. The ”[math]” point is the trivial vacuum, i.e. the zero-intensity configuration with only background lattice. The fluctuation equations about this point to quadratic order read

[TABLE]

where we have introduced the real variables $X_{1,3}=\Re\delta\Psi_{\pm},X_{2,4}=\Im\delta\Psi_{\pm},X_{5}=\delta E$ and

[TABLE]

The system (38) is degenerate at linear order thus we need a quadratic order expansion to analyze stability. The simplest approach is to construct a Lyapunov function for the equation (38). The function $V(X)=X^{2}$ is positive for and only for $X\neq 0$, and its derivative is

[TABLE]

which is strictly negative for $X$ nonzero if $f_{\pm}>0$ and $X_{5}>0$. However, we always have $X_{5}>0$ because $dX_{5}/dt$ in the full relaxation equations (37) has strictly negative right-hand side and $E$ grows monotonically from zero to $-I_{x}/(1+I_{x})$, and at any finite $t$ we have $E(t)-E(t=\infty)=X_{5}>0$. Thus the trivial equilibrium point is locally stable for $f_{+}>0,f_{-}>0$, i.e. $k>q^{2}$. It is much harder to construct the Lyapunov function for the global equations (37): in this case there are no additional symmetries and the stability of higher-dimensional systems is in general an extremely difficult topic. Thus there may well be regions far away from the ”[math]” point which do not flow toward it.

The ”$\pm$” pair is quite hard to study. All hope of expanding the system to second order and understanding the resulting complicated five-variable system is lost. This time, however, we can do a non-trivial first-order analysis as the system is non-degenerate and nicely reduces to the $(X_{1},X_{5})$ subsystem. Rescaling $X_{1}\mapsto(1+E_{0})^{-3/4}\left(I_{x}+E_{0}\left(1+I_{x}\right)\right)^{1/2}$ and $t\mapsto t\left(\left(1+E_{0}\right)/\left(I_{x}+E_{0}\left(1+I_{x}\right)\right)\right)^{1/4}$, the equation of motion for the ”$\pm$” point reads

[TABLE]

with $a_{\pm}$ being some (known) positive functions of $\Gamma,E_{0},I_{x}$ (independent of $k,q$). This is precisely the normal form for the Andronov-Hopf bifurcation arnold , and the bifurcation point lies at $k=-\Gamma E_{0}-q^{2}$. To remind, the bifurcation happens when the off-diagonal element in the linear term changes sign: the fixed point is stable when $a_{\pm}/(\Gamma E_{0}+k+q^{2})$ is positive. The sign of the nonlinear term determines the supercritical/subcritical nature of the bifurcation. Negative sign means the fixed point is stable everywhere before the bifurcation and is replaced by a stable limit cycle after the bifurcation (supercritical). Positive sign means the fixed point coexists with the stable limit cycle before the bifurcation and the $(X_{1},X_{5})$ plane is divided among their attraction regions; after the bifurcation there is no stable solution at all (subcritical).121212One should not take the stability in the whole $(X_{1},X_{5})$ plane in the supercritical case to seriously. We have expand the equations of motion in the vicinity of the fixed points and the expansion ceases to be valid far away from the origin.

In conclusion, stable ”$+$” equilibrium exists for $k>-\Gamma E_{0}-q^{2}$ where $E_{0}$ is best found numerically. Exactly the same condition holds for the ”$-$” point. For $k<-\Gamma E_{0}+q^{2}$, dynamics depends on the sign of the nonlinear term in (41): for the positive sign we expect periodically changing patterns. If the term is negative and the bifurcation is subcritical, various possibilities arise: the system may wander chaotically between the ”$+$” and the ”$-$” point, or it may end up in the attraction region of the ”[math]” point and fall onto the trivial solution with zero intensity. Naively, the attraction regions of the two fixed points (”$\pm$” and ”[math]”) are separated by the condition $-\Gamma E_{0}-q^{2}=q^{2}$, i.e. $q_{c}=\sqrt{-\Gamma E_{0}(\Gamma,\tau)/2}$, where we have emphasized that $E_{0}$ is in general non-universal. However, the actual boundary may be more complex, as our analysis is based on finite-order expansion around the fixed points, which is not valid far away from them.

The outcome is that the system generically has stable trivial and nontrivial (nonzero intensity) equilibria, in addition to time-dependent, periodic or aperiodic solutions. Numerical integration gives a similar picture of the stability diagram in Fig. 17. Numerically we find that the stability limit is $k>\Gamma-q^{2}$, i.e. $E_{0}\approx-1$. The region of applicability of our formalism lies in the top right corner of the diagram (nontrivial equilibrium), above $k\approx 1/L$. Formally, both $k$ and $q$ can be any real numbers. In practice, however, $k$ is discrete and its minimal value is of the order $1/L$. The spatial momentum $q$ lies between the inverse of the transverse length of the crystal (which is typically an order of magnitude smaller than $L$, i.e. minimal $q$ can be assumed equal to zero) and some typical small-scale cutoff which in our case is the vortex core size. We made no attempt to study the non-equilibrium behavior in detail or to delineate the boundary between the oscillatory and the chaotic regime since it is irrelevant for the main story of the paper.

From a practical viewpoint, the $\Gamma$-$I_{x}$ plane can be divided into two regions. One of them has a single stable ”$+$” or ”$-$” equilibrium or a $+\mapsto-$ limit cycle whose amplitude vanishes in the thermodynamic limit at all scales, i.e. for all $(k,\mathbf{q})$. This region can be legitimately described within the formalism of partition functions and equilibrium field theory. The second region flows toward the trivial fixed point and does not support vortices – this can also (trivially) be described by our formalism, as it always corresponds to the insulator regime, with no stable vortices. Thus the consistency check is that our method predicts no other phases in this region but insulator. In the third regime where long-term dynamics is either a limit cycle with amplitude of order unity, or chaos. This regime was studied in detail in some earlier publications (e.g. korn and references therein), and it cannot be reached within our present formalism.

B.2 Numerical checks

Now we complement the analytical considerations with numerical evidence that the phases described in the main text exist as long-term stable configurations. In Fig. 18 we show the time evolution of a vortex lattice in three different phases, where a visual inspection clearly suggests the system approaches equilibrium. In contrast, in Fig. 19 we see first a pattern that oscillates forever, i.e. follows a limit cycle (A), becomes incoherent (wandering chaotically over the unstable manifold, B), or dissipates away (reaching the [math]-fixed point), in the (C) panel. The loss of stability corresponds to an Andronov-Hopf bifurcation, as found earlier for non-vortex patterns in stabanal .

Dynamics can be most easily traced by looking at the numerically computed relaxation rate

[TABLE]

which is expected to reach zero for a generical relaxation process, where in the vicinity of an asymptotically stable fixed point $X\sim X_{\mathrm{eq}}+xe^{-rt}$, will be generically nonzero for a limit cycle or chaos, and will asymptote to a constant for the [math]-point, where $X_{\mathrm{eq}}=0$, so we get $(1/X)dX/dt\sim r$. Fig. 20 summarizes these possibilities. The black curves, corresponding to Fig. 18(a) and Fig. 18(c), show the situation which is in the focus of this work – the approach toward static equilibrium. The blue curve shows the limit cycle leading to periodic oscillations. The green curve corresponds to the chaotic regime with aperiodic dynamics and no relaxation, as in Fig. 19(b). Finally, the red curve corresponding to the pattern which radiates away in Fig. 19(b) reaches a constant value of $r$. In conclusion, the system shows roughly four classes of dynamics: fixed point, limit cycle, chaos and incoherence. Our work only covers the first of the four, but the bifurcation diagrams in the previous subsection give a good hint of the part of the parameter space which contains them, facilitating experimental or numerical verification.

Appendix C Perturbation theory and stability analysis

In this Appendix we develop the perturbation theory of the photorefractive beam system starting from the Lagrangian (4). The perturbation theory yields the criterion for the stability of the intensity patterns as they propagate along the $z$-axis. Formally, it is just the perturbative diagrammatic calculation of the propagator. This calculation explicitly excludes topologically nontrivial patterns and thus is somewhat peripheral for our main goal, understanding the vortex dynamics. But the general ways by which an envelope $\Psi_{\pm}$ can evolve along the $z$-axis and become unstable remain valid also for vortices. In particular, we will end up with a classification of geometrical symmetries of the intensity pattern $\Psi^{\dagger}\Psi$; the same symmetries are seen in vortex patterns and are an important guide for numerical and experimental work – how to recognize instabilities and also phases of the system.

Our system is strongly nonlinear, thus a naive perturbation theory about the ”trivial vacuum”, i.e. constant beam intensity is out of question. The right way is to perturb about a nontrivial solution, which approximates a stable pattern. This means we treat the light intensity as constant in ”time” $z$ but non-constant in space $(x,y)$. This is the hallmark of spatial dynamical solitons: they propagate with a constant profile along the $z$-axis and to a good approximation do not interact with each other and do not radiate korn . We thus write $\Psi=\Psi_{0}+\delta\Psi$, giving $\Psi^{\dagger}\Psi=I_{0}+F_{0}(\delta\Psi_{+}^{\dagger}+\delta\Psi_{+})+B_{0}(\delta\Psi_{-}^{\dagger}+\delta\Psi_{-})+\delta\Psi^{\dagger}\delta\Psi+O(|\delta\Psi|^{2})$ with $F_{0}^{2}+B_{0}^{2}=I_{0}$. The lowest-order Lagrangian for $\Psi_{0}$ now reads

[TABLE]

which determines the shape of the solution $\Psi_{0}(x,y)$ in the first approximation. The dynamical term with $\partial_{z}\Psi$ drops out (it is proportional to the equation of motion for $\Psi$). Nontrivial propagation in ”time” $z$ is obtained from second-order expansion of the potential which is given in the next Appendix in (52) and we will not copy it here. Varying the quadratic expansion with respect to the fluctuation $\delta\Psi$ gives the linearized equation of motion for $\delta\Psi$:

[TABLE]

where $\delta\Psi^{+}\equiv\delta\Psi^{\dagger},\delta\Psi^{-}\equiv\delta\Psi$. In homogenous ”spacetime” $(z,x,y)$ we can transform to momentum space in both transverse and longitudinal direction. In the transverse plane we get $(x,y)\mapsto(q_{x},q_{y})$ and $\Delta\mapsto-q^{2}$. The longitudinal coordinate or ”time” $z$ transforms as $z\mapsto k_{n}$ where $k_{n}=\pi n/L$, so the ”time” maps to discrete frequencies. The reason is of course that its domain is finite, corresponding to the crystal length $L$.

Now we can derive the bare propagator (Green’s function) of the fluctuating dynamical field $\delta\Psi$ by inserting the appropriate source $S(z)$ on the right-hand side of the equation (44). Normally, the source in the equation for the Green’s function is just the Dirac delta function but the counterpropagating nature of our beams imposes a ”two-sided” source:

[TABLE]

With this source (also Fourier-transformed in $z$), Eq. (44) gives the bare propagator $G^{(0)}_{\alpha\beta}$ for the fields $\delta\Psi^{\pm}_{\alpha\beta}$:

[TABLE]

The auxiliary matrices $A,B$ are defined as follows:

[TABLE]

where $P_{1}=(1/4)I_{0}\Gamma(1+\tau u+\tau E_{0})/(1+\tau u+I_{x}+I_{0})^{2},P_{0}=\Gamma-\Gamma(1+\tau u+\tau E_{0})/(1+\tau u+I_{x}+I_{0})$, and $S(k_{n})=\mathrm{diag}(1,e^{\imath k_{n}L})$.

Now we have the basic ingredient of the perturbation theory – the bare propagator. The self-energy correction $\Sigma$ of the propagator from the potential $V_{\mathrm{eff}}$ can be expanded in a power series over $\delta\Psi$, which gives an infinite tower of vertices. Simple combinatorial considerations give the expansion

[TABLE]

and the contraction over the internal indices of $\Psi^{\pm},\delta\Psi^{\pm}$ is understood. Now we can formulate the diagrammatic rules. We have two kinds of propagators, $G^{(0)}$ and its Hermitian conjugate. The mean-field values $\Psi_{0}^{\pm}$ are external sources. The term of order $(j_{1},j_{2},j_{3})$ contains $j_{1}+j_{3}$ propagator lines $G^{(0)}$ ($j_{1}$ of them ending with the source $\Psi_{0}$), and $j_{2}+j_{3}$ lines $\left(G^{(0)}\right)^{\dagger}$ ($j_{2}$ of them ending with a source $\Psi_{0}^{\dagger}$); altogether there are $j\equiv j_{1}+j_{2}+2j_{3}$ lines. The expansion has to be truncated at some $j$. Since the mass dimension of $\Psi$ is $1$, the $(j_{1},j_{2},j_{3})$-diagram has the scaling dimension $2-2(j_{1}+j_{2}+2j_{3})<0$, so all diagrams are irrelevant in the IR. This means we can make a truncation at small $j$.131313We do not worry about the UV divergences: we have an effective field theory and the UV cutoff is physical and finite. The leading terms are those where the order of the perturbation in $\delta\Psi^{\pm}$, which equals $j_{1}+j_{2}+2j_{3}$, is the smallest. This gives two classes of diagrams, one with $j_{1}=1,j_{2}=j_{3}=0$ and another with $j_{2}=0,j_{1}=j_{3}=0$. They contain a single external source and introduce the wavefunction renormalization, $G^{(0)}\mapsto ZG^{(0)}$, which does not influence the stability analysis. The four quadratic terms (with $(j_{1},j_{2},j_{3})=(2,0,0),(0,2,0),(1,1,0),(0,0,1)$) introduce a mass operator. Only the terms $(1,1,0)$ and $(0,0,2)$ are trivial (non-interacting); the other two are interacting as they contain $\left(\delta\Psi^{\pm}\right)^{2}$, and require the calculation of an internal loop, giving the dressed propagator:

[TABLE]

where the mass squared is a positive matrix, because the corresponding coefficients in (48) have positive sign (from the term $(-1)^{j_{1}+j_{2}+j_{3}}$ with $j_{1}+j_{2}+j_{3}=2$) and the integral of the bare propagator is also positive. Explicitly, it reads

[TABLE]

where the discrete ”frequency” $k_{n}$ is summed in steps of $\pi/L$. Other than the mass renormalization, the dressed propagator has the same structure as the bare one. Now we will consider what this means for the stability of the patterns.

C.1 The pole structure, stability and dispersion relations

Consider the poles of the propagator defined by the zeros of the eigenvalues of the matrix $G^{-1}_{\alpha\beta}(k_{n},q)$. The stable solution corresponds to the situation where the perturbation $\delta\Psi_{\pm}$ dies out along $z$, so the stability of the solution is determined by the condition that the pole in $q$ should have a non-positive imaginary part, i.e. that a small perturbation decays. The denominator depends on $k_{n},q$ solely through $k_{n}^{2},q^{2}$; it is linear in $k_{n}^{2}$ and quadratic in $q^{2}$. Therefore, each of the two eigenvalues $\lambda_{\pm}$ defines two pairs of opposite poles, $\pm q_{*+},\pm q_{*-},\pm q_{**+},\pm q_{**-}$. Out of these, two pairs are positive for all parameter values, so no imaginary part can arise, and we have either two pairs of centrally-symmetric imaginary poles, or one such pair, or none at all. We thus expect the sequence of symmetry-breaking transitions:

[TABLE]

Full circular symmetry is expected when there is no instability. With a single pair of unstable eigenvalues, we expect a square-like pattern with $\mathbb{C}_{4}$ symmetry, and with two pairs only a single reflection symmetry axis remains, yielding the group $\mathbb{C}_{2}$. Only in the presence of disorder in the background lattice intensity pattern $I_{x}$ can we expect the full breaking of the symmetry group down to unity, but this is an explicit breaking and is not captured by this analysis.

The dispersion relation for a typical choice of parameter values is represented in Fig. 21, where we plot the location of the pole $k(q)$ in the continuous approximation (interpolating between the $k_{n}$ values), with real parts of the pole in blue and imaginary in red. Since we have two pairs of opposite eigenvalues, the dispersion is $P$-symmetric in both $x,y$ and $z$ (remember that ”time” is really another spatial dimension), and any dispersion relation with a nonzero imaginary part will have a branch in the upper half-plane, i.e. an unstable branch. The only way out of instability is that the pole is purely real, i.e. infinitely sharp – this quasiparticle-like excitation signifies a solitonic solution. In Fig. 21, the dashed lines are drawn with the bare propagator $G^{(0)}$ and the full lines with the dressed propagator $G$, for the sets of parameter values. The perturbation always reduces the instability, i.e. the magnitude of the imaginary part of the poles – in (a, b) completely, resulting in zero imaginary part, and in (c, d) only partially. This reduction of instability likely explains the fact that linear stability analysis works extremely well for hypergaussian beams (which have most power at small values of $q$), as found in stabanal .

The fact that the imaginary region always lies at finite $q$ implies that the instability always starts at a finite scale, which corresponds to the behavior seen in the edge instability, which is shown, e.g., in Fig. 4. In order to understand the central instability, which starts from a single point, corresponding to $q\to\infty$, one needs to take into account also the higher order corrections from the potential (48) which, as we discussed, diverge at $q\to\infty$. While we always have a natural UV cutoff, it may happen that the corrections become large (though finite) before that UV scale is reached. We postpone a detailed account for the subsequent publication, and content ourselves to give only the diagram of the movement of the poles in the complex plane. Higher-order terms bring $q$-dependent corrections and break the inversion symmetry, resulting in the evolution of poles as in Fig. 22. The instability corresponds to the situations where at least one pole has a positive imaginary part, i.e. the first three situations in the figure. The last pattern, with no symmetry at all and two real poles, is stable (but not asymptotically stable, as there is no pole with non-zero negative imaginary part).

The analysis performed here is obviously incomplete, and we have contented ourselves merely to give a sketch of how the instabilities considered in the main text arise, as well as to formulate a perturbation scheme which allows one to study such phenomena. Further work along the lines of stabanal ; kivshar3 ; kivshar4 is possible by making use of our formalism, and we plan to address this topic in the future.

Appendix D Derivation of the vortex Hamiltonian from the microscopic Lagrangian

Starting from the vortex solution (6) we want to obtain an effective Hamiltonian for the vortex-vortex interaction. The task is to separate the kinetic term of the winding phase (with $\nabla\theta_{0\pm}=\sum_{i}Q_{\pm}\log|\mathbf{r}-\mathbf{r}_{i}|$) from (i) the intensity fluctuations $\delta\psi_{\pm}$ about some background value $\psi_{0\pm}$ and (ii) the non-vortex phase fluctuations ($\delta\theta_{\pm}$) in (6). The first task requires us to integrate out the amplitude fluctuations in the quadratic approximation. We first write $\Psi=\Psi_{0}+\delta\Psi$ and expand the Lagrangian to quadratic order:

[TABLE]

The zeroth-order (non-fluctuating) term $\mathcal{L}_{0}$ determines the intensity $I_{0}=\psi_{0}^{\dagger}\psi_{0}$, and produces the kinetic term for the vortex phase $\theta_{\alpha}$, which gives just two decoupled copies of the conventional XY vortex gas. The quadratic part $\mathcal{L}_{2}$ becomes quite involved when we separate the amplitude $\delta\psi$ and the phase $\delta\theta$. Inserting (6) into (52), one gets a quadratic action for $\delta\psi_{\alpha}$ and $\delta\theta_{\alpha}$. The rest gives a coupled quadratic action for the amplitude and phase fluctuations. Altogether, the Lagrangian is

[TABLE]

where $(\ldots)$ denote all terms of cubic or higher order in amplitude or phase fluctuations $\delta\psi_{\alpha},\delta\theta_{\alpha}$, and we have left out the constant terms independent of all field values. Primes denote the derivatives with respect to $r$. The first term defines the intensity fluctuations through $\psi_{\alpha}(r)$, and the second term (transformed through partial integration) yields the aforementioned conventional Coulomb gas of vortices after inserting the vortex solution from (6) for $\theta_{\alpha}$. The third term has the meaning of stiffness or mass matrix for intensity fluctuations and the last term gives rise to the coupling between the flavors, upon integrating out $\delta\psi$. The matrix $\hat{K}$ is

[TABLE]

The action is quadratic in $\delta\psi$; therefore, we know how to integrate it out and obtain an effective action depending only on phase fluctuations. To do that, we need to solve the eigenvalue equation for $\delta\psi$ obtained from (53), which reads

[TABLE]

and is solved by diagonalizing the system $\delta\psi_{\alpha}\mapsto\delta\chi_{\alpha}=U_{\alpha\beta}\delta\psi_{\beta}$ and reducing it to the Bessel equation:

[TABLE]

with $\nu=\sqrt{1/4+|\nabla\theta|^{2}}$, where $J_{\nu},Y_{\nu}$ are Bessel functions of the first and second kind, respectively. For well-defined behavior close to the vortex core (for $r\sim a$) we have $D_{\pm}=0$, and $C_{\pm}$ are arbitrary as the equation is linear. The eigenvalues $\lambda_{n\pm},n=0,1,ldots$ are obtained by requiring that the fluctuation decays to zero at the crystal edge $r=\Lambda$:

[TABLE]

where $j_{n}$ is the $n$-th zero of the Bessel function $J_{\nu}$. The values $\lambda_{n\pm}$ are impossible to express analytically in closed form, however it is not necessary for our purposes as $\Lambda\gg a$. The functional determinant obtained after integrating out $\chi_{\pm}$ is now expressed in terms of the eigenvalues:

[TABLE]

where $K^{\prime}_{\pm}=2K_{11}\mp K_{12}$. Now we are left with a solely phase-dependent quadratic Lagrangian:

[TABLE]

The final task is to integrate out the phase fluctuations, which is a trivial Gaussian integration, yielding:

[TABLE]

The resulting Lagrangian now depends only on the vortexing phases $\theta_{\alpha}$. The first term is carried from the original Lagrangian, and it does not mix the flavors. But the second term, stemming from the fluctuations, has nonzero mixed $\pm$ cross-terms. The quadratic derivative terms can be transformed by partial integration to the familiar Coulomb gas form of the XY model, with the same-flavor coupling which is already present in absence of fluctuations, and the coupling between the vortices of different flavors. Thus the existence of two beams together with the fact that amplitude and phase fluctuations do not decouple give us a richer system, with interaction between two vortex flavors. For future use it is more convenient to look at the vortex Hamiltonian $\mathcal{H}_{\mathrm{vort}}$ – the difference from the Lagrangian lies just in the sign of the term $V_{\mathrm{eff}}$. This finally yields the Hamiltonian (for Eq. (7), repeated here for convenience):

[TABLE]

with $r_{ij}\equiv|\mathbf{r}_{i}-\mathbf{r}_{j}|$, and the indices $1\leq i,j\leq N$ sum over all the vortices. The coupling constants $g,g^{\prime},g_{0},g_{1}$ are the result of integrating out the intensity fluctuations and in general are given by rather cumbersome (and not very illustrative) functions of $\Gamma,I_{0},\tau$. We give the expressions at leading order just for comparison with numerics:

[TABLE]

These expressions are used later to redraw the phase diagram in the space of physical parameters $\Gamma,I_{0},I_{x},L$.

Appendix E Multi-vortex mean-field theory

For a mean-field treatment of a system with multiple vortices, we start from the Hamiltonian (7) and introduce the order parameter fields in the following way. Denote the number of vortices with charge $(1,1)$ by $\rho_{2+}$ and the number of vortices $(1,-1)$ by $\rho_{2-}$; due to charge conservation, this means we also have $\rho_{2+}$ vortices of type $(-1,-1)$ and $\rho_{2-}$ vortices with charge $(-1,1)$. The number of single-charge vortices of type $(1,0)$ and $(0,1)$ is denoted by $\rho_{1+}$ and $\rho_{1-}$, respectively. Denote also $\rho_{2}\equiv\rho_{2+}+\rho_{2-}$ and $\delta\rho_{2}\equiv\rho_{2+}-\rho_{2-}$ (notice that $-\rho_{2}\leq\delta\rho_{2}\leq\rho_{2}$), and finally $\rho_{1}\equiv\rho_{1+}+\rho_{1-}$. We insert this into the vortex Hamiltonian $\mathcal{H}_{\mathrm{vort}}$, and assume that the long-ranged logarithmic interaction $\log r_{ij}$ justifies the mean-field approximation: for $i\neq j$ we can approximate $\log r_{ij}\sim\log\Lambda$, assuming that average inter-vortex distance is of the same order of magnitude as the system size. For the core energy, we know that $g_{0},g_{1}\sim\log(a/\epsilon)\sim-\log\epsilon\sim\log\Lambda$, where in the last equality we have assumed that the UV cutoff $\epsilon$ is of similar order of magnitude as the inverse of the IR cutoff $1/\Lambda$, which is natural.141414Nevertheless, this is clearly not a rigorous argument. Our mean-field calculation is somewhat sketchy and merely assumes that the long-range interactions can safely be modeled as a uniform vortex charge field. Thus all terms are proportional to $L\log\Lambda$ and we can write

[TABLE]

We use the notation $\beta\equiv L$, to emphasize the analogy with the free energy of spin vortices, where $\beta$ is the inverse temperature. The analogy is purely formal as our system is not subject to thermal noise. Now the ground state is determined by minimizing the free energy, i.e. the effective action of the system. Notice that $\mathcal{F}_{\mathrm{mf}}$ is linear in the fields $\rho_{2},\delta\rho_{2},\rho_{1}$ so the optimal configurations have either $\mathcal{F}_{\mathrm{mf}}=0$ or $\mathcal{F}_{\mathrm{mf}}\to-\infty$, and the mean-field densities $\rho_{1,2}$ are either zero, or arbitrary (formally infinite). This is a well-known property of the 2D Coulomb gas and has to do with the fact that (assuming the cutoff dependence has been eliminated) this system is conformal invariant in the insulator phase, so all finite densities $\rho$ are equivalent: there is no other scale to compare $\rho$ to. Likewise, the prefactor $\beta$ can be absorbed into the definition of the coupling constants $g,g^{\prime}$ and thus is not an independent parameter (this is well-known also from the single-flavor case). Minimizing (63) is an elementary exercise and we find again four regimes, corresponding to the four phases we guessed based on the single-vortex free energy $\mathcal{F}^{(1)}$:

For $A>0,A>|B|$ the minimum is reached for $\rho_{2}=\delta\rho_{2}=\rho_{1}=0$. In the ground state there are no vortices at all – the system is a vortex insulator. 2. 2.

For $B>0,A+B<0$, giving $g+g^{\prime}<1,g^{\prime}>0$, the free energy has its minimum for $\rho_{2}>0$ and $\delta\rho_{2}=-\rho_{2}$ (notice that $-\rho_{2}\leq\delta\rho_{2}\leq\rho_{2}$). This means $\rho_{2+}=0,\rho_{2-}>0$, so opposite-charged vortices $(Q,-Q)$ proliferate, and the system is dominated by the interactions between the charges. This is the frustrated vortex insulator regime. Since $g^{\prime}<0$, the single-charge vortices (density $\rho_{1}$) are suppressed. 3. 3.

For $B<0,A+B<0$, i.e. $g+g^{\prime}<1,g^{\prime}<0$ the minimum is reached for $\rho_{2}=\delta\rho_{2}>0$, i.e. $\rho_{2-}=0$ so the vortices $(Q,Q)$ can proliferate. However, since $g^{\prime}<0$ there is also nonzero single-flavor density $\rho_{1}$ and the proliferation of vortices $(Q,0)$ and $(0,Q)$ which generically dominate over two-flavor vortices. This is the conductor phase, with mostly single-flavor vortices (as in the standard XY model). 4. 4.

The point $A=B=0$ is special: naively, from (63) arbitrary nonzero $\rho_{1},\rho_{2},\delta\rho_{2}$ are allowed. Of course, higher-order corrections will change this but the energy cost of vortex formation will generically be smaller than in previous phases. This is the vortex perfect conductor phase. In the mean-field approach it looks like a single point, but that will turn out to be an artifact of the mean-field approach: for small nonzero $A,B$ the system still remains in this phase.

In terms of the original parameters $g,g^{\prime}$, one sees the insulator phase is given by $g+g^{\prime}>1$, and the conductor and the FI are separated by the line $g^{\prime}=0$. We can now sketch the phase diagram, which is given in Fig. 2(a), side by side with the more rigorous diagram obtained by the RG flow, in the subsubsection IIIB2.

Appendix F Counter-propagating boundary conditions

In the derivation of the vortex Hamiltonian and its RG analysis, we have pulled under the rug the treatment of the CP boundary conditions: the effective Hamiltonian (and consequently the partition function and the phase diagram) depends solely on the bulk configuration, and nowhere can one see the fact that $\Psi_{+}(z=0;\mathbf{r},t)$ and $\Psi_{-}(z=L;\mathbf{r},t)$ are fixed. Now we will explicitly show that these boundary conditions are irrelevant in the RG sense, i.e. they contribute additional, boundary terms to the effective Hamiltonian, but these terms do not change the fixed points to which the solution flows.

The full Hamiltonian with correct CP boundary conditions is obtained by adding the $F$-source at $z=0$ and the $B$-source at $z=L$ to the Lagrangian $\mathcal{L}$ from (4) or, equivalently, to the equations of motion. The sources impose the conditions $F(z=0;x,y;t)=F_{0}(x,y)$ and $B(z=L;x,y;t)=B_{0}(x,y)$ so they equal

[TABLE]

and the full Lagrangian is

[TABLE]

Unlike the Dirac delta source (45) for the Green’s function, now the source has nontrivial dependence on transverse coordinates. Now we can insert the vortex solution (6) in both $\Psi_{\pm}$ and $F_{0},B_{0}$ and repeat the steps from the subsequent derivation. The vortex charges in $F_{0}$ can be denoted by $\vec{P}_{i^{\prime}}^{(+)}\equiv(P_{i^{\prime}+},0)$ and $\vec{P}_{i^{\prime}}^{(-)}\equiv(0,P_{i^{\prime}-})$; by definition, the $+$ component of $B_{0}$ as well as the $-$ component of $F_{0}$ are zero and thus carry no vorticity. The primed indices refer to the vortices in the input beams, and the non-primed, like before, to the bulk vortices. Notice the source term changes sign upon performing the Legendre transform, appearing as $-J_{+}\Psi_{+}-J_{-}\Psi_{-}$ in the Hamiltonian.

Now we will check if the RG flow of the Hamiltonian with boundary terms is affected by the sources. In the notation introduced above, the total vortex Hamiltonian is

[TABLE]

Notice there is no source-source interaction: same-flavor interaction cannot exist as $P_{i^{\prime}-}=P_{j^{\prime}+}=0$, and cross-flavor interaction does not exist as $J_{+}$ and $J_{-}$ exist at different $z$ values, i.e. the cross-term would be proportional to $\delta(z)\delta(z-L)$ and thus vanishes. The presence of sources breaks the spatial homogeneity, complicating the traces (integrals over the positions of virtual vortex-antivortex pairs), but does not change the main line of the calculation. The fluctuation of the partition function due to vortex pair creation is now

[TABLE]

We have denoted $C\left(\vec{Q}_{1},\mathbf{R}_{1};\vec{Q}_{2},\mathbf{R}\right)\equiv\left(g\vec{Q}_{1}\cdot\vec{Q}_{2}+g^{\prime}\vec{Q}_{1}\times\vec{Q}_{2}\right)\log R_{12}$, and the coupling to the sources is encapsulated in the function

[TABLE]

with $z_{+}=0,z_{-}=L$. The coupling constants $g_{\pm},g^{\prime}_{\pm}$ are obtained from $g,g^{\prime}$ in (67) by replacing

[TABLE]

with $I_{+}=|F_{0}|,I_{-}=|B_{0}|$. Now the exponential of the extra term with sources also needs to be expanded in $r_{12}$, as the combinations $2\mathbf{r}=\mathbf{r}_{1}+\mathbf{r}_{2},\mathbf{r}_{12}=\mathbf{r}_{1}-\mathbf{r}_{2}$ do not decouple anymore and exact integration is impossible. After writing $D_{j^{\prime}}^{\pm}(\vec{q},\mathbf{r}_{1,2})=D_{j^{\prime}}^{\pm}\left(\vec{q},\mathbf{r}\right)\pm\mathbf{r}_{12}\cdot\nabla D_{j^{\prime}}^{\pm}\left(\vec{q},\mathbf{r}\right)+\ldots$ and similarly for $C$, we notice first that the zeroth-order terms from $D_{j^{\prime}\pm}$ cancel out: $D_{j^{\prime}}^{\pm}\left(\vec{q},\mathbf{r}\right)+D_{j^{\prime}}^{\pm}\left(-\vec{q},\mathbf{r}\right)=0$. Then we get (at quadratic order in $r_{12}$ in the integrand):

[TABLE]

The integral $I_{mn}$ is the term with the contribution of order $r_{12}^{m}$ from the second line in the integrand and with the contribution of order $r_{12}^{n}$ from the third line. The integrands in $I_{mn}$ are thus of order $m+n$ in $\mathbf{r}_{12}$, $m$ coming from the expansion of $D_{j^{\prime}}^{\pm}$ and $n$ from the expansion of $C$. By homogeneity $I_{01}=0$ and $I_{02}$ is the same integral that appears in absence of sources, whose calculation was used in obtaining (13) and which gives the right-hand side of the RG flow (14). The remaining integral $I_{11}$ is the new ingredient, and the only one which depends on the sources. Representing it as

[TABLE]

we compute the integral $\tilde{I}_{j^{\prime}\alpha}$ in polar coordinates:

[TABLE]

where $\theta_{j^{\prime}\alpha}$ is the angle between $\mathbf{r}_{j^{\prime}\alpha}$ and $\mathbf{r}_{12}$. Assuming the RG scale changes as $\Lambda_{1}=\Lambda,\Lambda_{2}=\Lambda(1+\ell)$, for small $\ell$ we can expand the integrand getting:

[TABLE]

The complicated dependence on the positions of the sources disappears completely in the first order in $\ell$!151515The additional assumption is that $\Lambda>r_{j^{\prime}\alpha}$ so the integrand contains no poles; this is clearly justified as $\Lambda$ is the length cutoff. Altogether, comparing the outcome of (70) to the original Hamiltonian (66), we see that the renormalization of the bulk interaction between $\vec{Q}_{1}$ and $\vec{Q}_{2}$ is unaffacted by the sources, given as before by the $I_{02}$ term, and the source-bulk coupling renormalizes with a strictly negative shift (as $\tilde{I}_{j^{\prime}\alpha}=2\pi>0$). The flow equations for $g,g^{\prime}$ couplings are unchanged, being the same as in (14). The bulk-to-boundary couplings $g_{\pm},g^{\prime}_{\pm}$ have the flow equations

[TABLE]

where $N=\sum_{j^{\prime}}\sum_{\alpha}1$ is the total vorticity of the sources. This obviously flows to $g_{\pm},g^{\prime}_{\pm}=0$.

Intuitively, one may wonder how come that such an important thing as the CP geometry has no bearing on the vortex dynamics; surely the behavior of a co-propagating system would be expected to differ from a counter-propagating system. The answer is that the CP geometry does enter our calculations – the $B$-beam has an extra minus sign in the equations of motion (1) (alternatively, in the Lagrangian in Eq. 4); equivalently, the symmetry group of the effective potential in the Lagrangian is $SU(1,1)$, not $SU(2)$ as it would be for two co-propagating beams. Finally, let us emphasize again that in the numerical simulations we directly solve the propagation equations (1) together with (2), i.e. we directly take into account the CP boundary conditions – no analytical approximations whatsoever are used in the numerics, and no use is made of the effective vortex Hamiltonian.

Appendix G Order parameters and RG analysis of the CP vortices in the presence of disorder

G.1 Saddle-point solutions

We start by rewriting the replicated partition function $\bar{\mathcal{Z}}^{n}$ in terms of $p_{\alpha},q_{\alpha\beta}$ and inserting the constraints which encapsulate their definition in Eq. (19):

[TABLE]

We have five constraints, $\lambda_{(\mu\nu)}^{++},\lambda_{(\mu\nu)}^{--},\lambda_{(\mu\nu)}^{+-}=\lambda_{(\mu\nu)}^{-+},\lambda_{(\mu)}^{+},\lambda_{(\mu)}^{-}$, for the corresponding five order parameters in (19). We can denote

[TABLE]

We will also sometimes leave out the replica indices $\mu,\nu$ to avoid cramming the notation too much. Now we can first integrate out the vortex degrees of freedom $Q_{i\alpha}^{(\mu)}$ from (18) to get the effective action

[TABLE]

The saddle-point equations for the constraints give the constraints in terms of the expectation values $q^{(\mu\nu)},p^{(\mu)}$. Luckily, the equation for $\vec{\lambda}$ is easy:

[TABLE]

so we immediately solve $\vec{\lambda}=\hat{K}\vec{p}$. Now plugging this into the equations for the three remaining constraints yields

[TABLE]

We have denoted $X=\det\lambda_{++}=\det\lambda_{--},Y=\det\lambda_{+-}$ (these have a well-defined limit for $n\to 0$). It is trivial to write $\lambda_{\pm\pm},\lambda_{+-}$ from the above expressions, and we can feed the solutions for all the constraints into the effective action and then solve the saddle-point equations for the order parameters $p_{\pm},q_{++},q_{--},q_{+-}$. Full equations are too complex to be solved, even approximately. We will simplify the problem with the following reasoning. The sums over single-replica order parameters generically scale as $\sum_{\mu}p^{(\mu)}_{\pm}\sim\sum_{\mu}\left(p^{(\mu)}_{\pm}\right)^{2}\sim n$, whereas the double-replica parameters have $\sum_{\mu,\nu}q^{(\mu\nu)}_{\alpha\beta}\sim n^{2}$. This means that in the limit $n\to 0$ the $p_{\pm}$-terms dominate over $q_{\alpha\beta}$-terms. Therefore, if $p_{\pm}\neq 0$ we can disregard the quantities $q_{\alpha\beta}$ or expand in a series over them, simplifying the equations significantly. Only if the replica symmetry breaking imposes $p_{\pm}=0$ (not every replica-symmetry-breaking configuration does so), the $q_{\alpha\beta}$ order parameters are significant, and the saddle-point equations with $p_{\pm}=0$ are again approachable.

Consider first the case $p_{\pm}=0$. After some algebra, the effective action is now

[TABLE]

Consider first the ansatz when the $q^{\pm\pm}$ fields are nonzero, whereas the mixed-flavor field $q^{+-}$ is zero. In this case, the second term in (82), coming from the determinant $\hat{K}$, simplifies further and we get the saddle-point equation

[TABLE]

which is the same as for the infinite-range spin-glass Ising model cardy ; spinglass . One obvious solution is $q^{\pm\pm}=q^{+-}=0$, the completely disordered system with no vortex proliferation – the familiar insulator phase. It is easy to check that this is indeed a minimum of the effective action $S_{\mathrm{eff}}$. There is also a replica-symmetric but nontrivial solution

[TABLE]

which yields the solution $\mathcal{Q}_{0}^{\pm\pm}=1-1/(\beta\sigma^{\pm\pm})$. However, this solution is unstable and is not observable. A stable non-trivial solution is obtained if the replica symmetry is broken. The ansatz is well-known from the spin-glass literature (e. g., spinglass ) and has a $\rho\times\rho$ matrix $\hat{Q}^{\pm\pm}$ on the block-diagonal and the constant zero elsewhere, with

[TABLE]

Eq. (83) suggests that $\mathcal{Q}_{1}^{\pm\pm}>0$ for sufficiently large $\beta$, i.e. small $L$. However, no analytical solution for the elements $\mathcal{Q}_{1}^{\pm\pm}$ exists and they have to be solved for numerically, by plugging in the solution into the effective action and minimizing it. This is an easy task (for chosen values of the parameters and disorder statistics) but we will not do it here as we do not aim at quantitative accuracy anyway; we merely want to sketch the phase diagram. Now if the third field $q^{+-}$ is nonzero, it satisfies the same equation as (83) just with $\sigma_{++}^{2}\mapsto 2\sigma_{+-}^{2}$. The three combinations of nonzero order parameters correspond to the three familiar phases – $q^{\pm\pm}\neq 0$ is the conductor, $q^{+-}\neq 0$ is the frustrated insulator and $q^{\pm\pm},q^{+-}\neq 0$ is the perfect conductor.

The solutions with $\vec{p}\neq 0$ yield new physics. In this case, we have at leading order $\lambda_{\pm\pm}=-2X/(X^{2}-Y^{2})(p^{\pm})^{-2},\lambda_{+-}=-2Y/(X^{2}-Y^{2})(p^{+}p^{-})^{-1}$, so the effective action is

[TABLE]

giving the saddle-point equations

[TABLE]

which easily gives

[TABLE]

with $s_{1,2}\in\{\pm 1\}$. The solution is the same for every $\mu$ and $p^{\pm}_{(\mu)}=\left(p^{\pm},p^{\pm},\ldots,p^{\pm}\right)$. Now depending on the sign of the determinant $J_{0}^{+}J_{0}^{-}-\left(J_{0}^{+-}\right)^{2}$ the solutions for different $s_{1,2}$ may be minima or saddle points. In either case, we have a phase with nonzero local charge density, which is the meaning of $\vec{p}$. If there are multiple minima we call this phase vortex glass. The reader may argue that true glass should satisfy more stringent conditions and that our phase is not a true glass. Depending on the viewpoint this may well be accepted and we use the term ”glass phase” merely as shorter and more convenient than ”phase with power-law correlation decay, no long-range order and frustrated free energy landscape”. The phase with a single minimum will be called charge density wave, as it has a unique ground state configuration yielding macroscopically nonzero charge density, i.e. it has a true long-range order. On the other had, with multiple minima the replica-averaged charge density sums to zero. The landscape, i.e. the effective action of the system for given replica ($\mu$) as a function of $p^{\pm}_{(\mu)}$, is given in Fig. 23 as the density map of the function $S_{\mathrm{eff}}(p^{+},p^{-})$ dependence for $J_{0}^{+}=-J_{0}^{-}=1$ (glass phase, A) and $J_{0}^{+}=J_{0}^{-}=1$ (charge density wave, B). We see that the glassy phase shows two inequivalent minima in each replica, with $s_{1}=-s_{2}=\pm 1$ in Eq. (88), so the total action, the sum of actions of all replica subsystems, can have one and the same value for many configurations, the definition of a highly frustrated system, one of the reasons we dub this phase glass. The charge density wave only has a single minimum for $s_{1}=s_{2}=1$.

G.2 RG flow equations

The starting Hamiltonian is the same as in (18). Now we will write it out more explicitly, keeping the distance-dependent parts:

[TABLE]

We have denoted the elements of $J_{0}$ by $J_{0}^{++}=J_{0}^{--}=\bar{g}_{c},J_{0}^{+-}=J_{0}^{-+}=\bar{g}^{\prime}_{c}$ (the bars over the letter remind us that these are disorder-averaged values). The fluctuation of the partition function is completely analogous to the clean case, only it has the additional non-local quartic term. It can again be expanded over $r_{12}$ as in (12) but the quartic term contains no small parameter for the power series expansion and has to be kept in the exponential form. Starting from the expression for the fluctuation analogous to the clean case (12), we get

[TABLE]

We have used the notation:

[TABLE]

Now we trace out the fluctuations first by integrating over $r$ and doing some simple algebra:

[TABLE]

The next step is the summation over all possible $\pm 1$ charges of virtual vortices $\vec{q}^{(\mu)},\vec{q}^{(\nu)}$ (the two replica indices mean two summations from $1$ to $n$), which requires quite some algebra. The renormalized partition function $\bar{\mathcal{Z}}^{n}$ finally gives the RG flow equations:

[TABLE]

As discussed in the main text, the fixed point must lie either at $y=0$ or $y\to\infty$, depending on the magnitude of $g+g^{\prime}+\beta^{2}\sigma^{2}$. For $y\to 0$, three clean fixed points remain, which flow to zero disorder: these correspond to PC, FI and conductor. The disordered fixed point also has $y\to 0$ but the disorder is nonzero: this is the CDW phase from the mean-field analysis, the dirty analog of the insulator. Finally, when $y\to\infty$ and nonzero $\sigma^{2}$ at the fixed point we expect glassy behavior.