Dirac equation on a square waveguide lattice with site-dependent coupling strengths and the gravitational Aharonov-Bohm effect
Christian Koke Changsuk Noh Dimitris G. Angelakis

	
	
	
	
	christian.koke@campus.lmu.de
	Arnold Sommerfeld Center, Ludwig Maximilian University of Munich, Munich 80333, Germany
	
	
	cnoh@knu.ac.kr
	Department of Physics, Kyungpook National University,  Daegu 41566, Korea
	
	
	
	Centre for Quantum Technologies, National University of Singapore, Singapore 117542
	School of Electrical and Computer Engineering, Technical University of Crete, Chania 73100, Greece
	
	
	
		The main objective of this work is to present a theoretical proposal for an implementation of the (2 + 1)-dimensional Dirac equation in classical gravitational and electromagnetic backgrounds in a two-dimensional	waveguide array. For this, a framework for achieving site-dependent effective coupling constants in two-dimensional
		waveguide arrays is developed. Implementability of the Dirac equation under the proposed scheme puts minor demands on gauge and spacetime backgrounds; however, a wide array of physical spacetimes, such
		as all vacuum and static solutions, prove to be implementable. As an interesting and instructive example, we discuss a tabletop realization of the gravitational Aharonov-Bohm effect: After devising a thought experiment in which signatures of the gravitational Aharonov-Bohm effect could be detected, we briefly discuss how the analogue of such a setting can in principle be implemented using the proposed waveguide setup.
	
	
	
	Dirac equation on a square waveguide lattice with site-dependent coupling strengths and the gravitational Aharonov-Bohm effect
    Dimitris G. Angelakis
    Received May 29, 2019; revised July 5, 2019; accepted
July 10, 2019. Date of publication July 31, 2019 
 date of current version
August 21, 2019
====================================================================================================================================================



§ INTRODUCTION

	Physical effects stemming from global geometric or topological properties of a system have long been investigated in physics. A premier example of such an effect is the celebrated Aharonov-Bohm effect in  quantum mechanics <cit.>. 
	Similar phenomena can occur in 2+1 dimensional gravity, where localised sources do not influence the curvature of the manifold beyond their point of localisation, but may effect the global geometry <cit.>.
	In this work we propose a classical optical simulation of a gravitational analogue of the Aharonov-Bohm effect for Fermions in 2+1 dimensions <cit.>.
	
	There exists a multitude of  proposals and experimental demonstrations of interesting physics in coupled waveguide arrays
	<cit.>.
	Optical simulations of the (1+1)-dimensional Dirac equation in binary waveguide arrays have long been proposed <cit.> and experimentally demonstrated <cit.>.
	Simulations of Dirac dynamics in curved spacetime have also been investigated. 
	For approaches in graphene see Refs. <cit.>; investigations using cold atoms, e.g. Refs. <cit.>, are also conducted.
	For the general simulation of the (1+1)-dimensional Dirac equation in curved spacetime in coupled waveguide arrays, see <cit.>.
	We now expand on this work by elaborating on how the (2+1)-dimensional Dirac equation may be simulated in coupled waveguide arrays and then show how one can incorporate classical gravitational and electromagnetic background fields.
	
	This article is organized as follows: After discussing the types of metrics and gauge connections we will consider, we recall the basics of the (2+1)-dimensional Dirac equation coupled to background fields in section <ref>. Subsequently, generalizing an idea proposed in Ref. <cit.>, we show how one can achieve site dependent effective coupling constants in two-dimensional coupled waveguide arrays in section <ref>. Section <ref> then explains how one can use this to implement the Dirac equation in square lattice waveguide arrays. Finally, in section <ref> we develop a scenario in which the 'phase shift' of the gravitational Aharonov-Bohm effect can be observed and discuss its potential realization in a coupled waveguide array.	
	

§ THE MINIMALLY COUPLED DIRAC EQUATION IN CURVED SPACETIME 

	In what follows, we shall always assume a global manifold chart.
	

 §.§ General considerations

	Our starting point is the coordinate expression of the minimally coupled Dirac equation in curved spacetime:
	
    [ iγ^μ∇_μ -m]ψ(x) = 0

	Here we are using the vielbein e^μ_a(x),
	to transform the local gamma matrices γ^a according to
	
    γ^μ(x) = e^μ_a(x)γ^a.

	The covariant derivative is given by
	
    ∇_ν = ∂_ν + Ω_ν + i A_ν,

	with a  gauge connection  A_ν and where  Ω_ν arises from the spin connection <cit.>.
	For a brief introduction to the Dirac equation in curved spacetime, the reader might want to consider the Appendix.
	
	
	
	For reasons that will become apparent shortly, we assume that the metric is of the form
	
    ds^2= e^2Φ(x,y,t)dt^2-e^2Ψ(x,y,t)(dx^2+dy^2),

	and that 
	
    Φ(x,y,t) - Ψ(x,y,t) = ρ(x,y),

	where ρ is independent of time.
	These restrictions still allow for all vacuum solutions, since, as argued in Appendix, vacuum solutions can be taken to be conformally flat and hence we may even choose Φ = Ψ.  Furthermore we note that any static spacetime may be represented in this form <cit.>. 
	For the gauge connection, we demand that the vector part may be gauged away, i.e. that we may write
	
    A_μd^μx=   ϕ(x,y,t) dt + 
       ∂_xΛ(x,y,t) dx + ∂_y Λ(x,y,t) dy.

	
	
	
	Under these assumptions, after choosing γ^0 = σ_z, γ^1=  - iσ_x and γ^2=  iσ_y, 
	and first neglecting the gauge part, the Dirac equation can be written as 
	
    i(∂_t+Ψ_t)ψ= -iσ_ye^Φ-Ψ(∂_x + 1/2Ψ_x +1/2Φ_x )ψ
       -iσ_xe^Φ-Ψ(∂_y + 1/2Ψ_y +1/2Φ_y )ψ  +σ_ze^Φ m ψ.

	See Appendix <ref> for calculations leading up to this.
	Recalling that Φ-Ψ = ρ and reinstating the gauge connection yields 
	
    i(∂_t+Ψ_t  +iϕ)ψ   = +σ_ze^Φ m ψ
    
    			-iσ_ye^ρ   (∂_x + 1/2Ψ_x +1/2Φ_x +i∂_x Λ)ψ
    
    			-iσ_xe^ρ   (∂_y + 1/2Ψ_y +1/2Φ_y +i∂_y Λ)ψ.

	Upon defining
	
    χ := e^[1/2(Φ + Ψ)+iΛ]ψ ,

	equation (<ref>) reduces to 
	
	
	
	
	
	
    ∂_tχ=   -ie^ρ(σ_y∂_x +σ_x∂_y )χ
       +(σ_ze^Φ m + 1[ϕ- ∂_t Λ] ) χ.
 
	
	We have assumed the metric to be of the form (<ref>) because this is precisely the form for which the resulting Dirac equation is implementable with our proposed scheme in a two-dimensional coupled waveguide array, below. As will be explained in section <ref>, our scheme allows implementation of partial differential equations precisely of the form
	
    ∂_tχ=   -ie^ρ(x,y)(σ_y∂_x +σ_x ∂_y )χ
       +(σ_z ĥ(x,y,t)  + 1h̃(x,y,t) ) χ,

	
	in a two dimensional square photonic lattice.
	
	To see how (<ref>) constrains metric and vector potential, it is expedient to consider the gauge and gravity cases separately.
	For the gauge part, note that while σ_z occurs independently of any differential operator in equation (<ref>), the Pauli matrices σ_x,y are always multiplied by partial derivatives. Thus we see that we have to be able to gauge away the spatial  part of the vector potential. If this were not the case there would be terms in which the Pauli matrices σ_x,y would appear independently of partial derivatives, since the gauge covariant derivative is given by ∇_μ = ∂_μ + i A_μ,
	
	
	
	which is contracted with the gamma-matrices {γ^μ}. 
	To be able to carry out a field redefinition like (<ref>) that cancels the x and y components of the gauge potential in the differential equation, we indeed have to be able to express the spatial part as ∂_xΛ dx + ∂_y Λ dy for some function Λ.
	
	For the gravity part, note that since each derivative term is multiplied by a single Pauli matrix and not a linear combination of them (at fixed position and time), the metric must necessarily be diagonal for our choice of gamma matrices. If it were not, the off-diagonal elements of the vielbein would give rise to terms such as  γ^a ∂_b,  (a≠ b).
	To ensure that possible terms arising from the gravity part, where σ_x,y Pauli matrices would appear independently of derivatives, can be absorbed into a field redefinition such as (<ref>), we have to demand that the two nonvanishing entries in the spatial sector of the metric are the same.
	
	Note that this then also implies that we may without loss of generality consider the functions B and C in (<ref>) to be identically zero.
	Lastly, we note that after bringing the Dirac equation into its Hamiltonian form, demanding that the function ρ= ρ(x,y) that modifies the strength of the derivative terms be time-independent is equivalent to assuming  Φ-Ψ = ρ to be independent of time. This is precisely the demand in (<ref>).
	
	
	

 §.§ Example: Conical space and the gravitational Aharonov-Bohm effect

	Before ending this section let us consider a specific example. Setting the gauge connection to zero, we assume that the metric is given by
	
    ds^2    = dt^2 - r^-2 Δ(dr^2+r^2 dθ^2)
       =dt^2 - (√(x^2+y^2))^-2 Δ(dx^2+dy^2).
 
	A short discussion on the underlying geometry is provided in Appendix  <ref>. 
	Defining
	
    f:=1/g:= e^-Ψ = (√(x^2+y^2))^Δ,

	the Dirac equation in this spacetime can be written as 
	
    ig∂_t ψ =    - 1/√(g)iσ_y(∂_x + g_x/2g)ψ -1/√(g)iσ_x(∂_y + g_y/2g)ψ   
       + g σ_z ψ.

	Upon redefining χ:= √(g)ψ this further reduces to
	
    i ∂_t χ =if( σ_y ∂_x + σ_x ∂_y) χ+ mσ_z χ.

	
	Let us now explore what happens to a spinor that is parallel transported along a curve in this spacetime. Let us assume that the curve is parametrized by ξ, in which case the tangent vector to the curve be given by ξ̇.
	To find the expression for a parallel transported spinor, we have to solve the equation
	
    ∇_ξ̇ψ = ξ̇^μ∇_μψ = 0.

	With the help of (<ref>) and by writing ψ =Vψ_0, where V:=V_ξ(τ) is defined to be the operator that parallel transports a spinor along the curve ξ, we may write this as 
	
    -dV/d τ = i σ_z/2(ξ̇^y Ψ_x - ξ̇^x 	 Ψ_y) V,

	where ξ̇^x/y is the x/y component of ξ̇. 
	
	Next we note that in the case at hand we have Ψ(t,x,y) = Ψ(r), with r= √(x^2+y^2).  Hence
	
    dV/d τ V^-1 = i ∂_r Ψσ_z/2r (ξ^xξ̇^y - ξ^yξ̇^x),

	with r^2 = (ξ^y)^2 + (ξ^x)^2.
	
	Assuming that we parallel transport around a concentric circle, we may write  ξ^1= rcos(τ), ξ^2= rsin(τ), and hence write the above equation as 
	
    dV/d τ V^-1 = i r∂_r Ψσ_z/2 (-sin^2(τ) -cos^2(τ)).

	Noting that
	
    r ∂_r Ψ = -Δ,

	one readily obtains 
	
    V(τ) = exp( Δiσ_z/2τ).

	It can in fact be shown that the assumption that we are on a concentric circle is not necessary and this expression is valid for general curves <cit.>.
	We thus note that parallel transport of a spinor around the center once, amounts to a multiplication with
	
    U(2 π)= cos(Δπ) + i σ_zsin(Δπ).


In section <ref>, we will propose an experiment to observe this non-trivial rotation in a tabletop experiment. 



§ EFFECTIVE COUPLING CONSTANTS IN COUPLED WAVEGUIDE ARRAYS

We want to construct an implementation of (<ref>) in a coupled waveguide array. 
Our starting point for this endeavour is the usual description of the array in terms of coupled mode equations.
For a 2D square array these are given by, 

    i∂c̃_ab/∂ z =    k_0 ( c̃_(a+1)b + c̃_(a-1)b +c̃_a(b-1) +c̃_a(b+1))
        + σ_ab(z) c̃_ab.

assuming only nearest neighbour and negligible diagonal couplings  <cit.>. 
Here z is the Cartesian coordinate along the waveguides,   c̃_ab denotes the amplitude in the waveguide at site (ab), σ_ab(z) is the difference between the refractive index of the waveguide and the bulk refractive index at site (ab) and position z.
The physical coupling strength determined by the overlap between the transverse components of the modes in adjacent waveguides is denoted by k_0.
Neglecting diagonal couplings is justified if the overlap of the evanescent tails between waveguides that are separated by √(2)δ_s is negligible compared to the overlap for wave guides that are separated by δ_s. In Appendix <ref> we derive a scheme to effectively cancel any diagonal couplings, in a general setting. 

To allow for site-dependent coupling constants, we generalize the idea proposed in <cit.> and introduce refractive index modulations at every site:

    σ_ab(z) = ∑_i=1^ΛA_ab^i/2 cos(α_i z)

Then we apply the transformation 

    c̃_ab = c_abexp(i  ∑_j=1^Λ A^j_ab/2 α^jsin (α_j z) ),

which cancels the sinusoidal terms in the coupled mode equations and obtain

    i∂ c_ab/∂ z =       k^*_[a↔ (a+1)]b c_(a+1)b 
    	+	  k_[(a-1)↔ a]b c_(a-1)b
    
    	+	     k^*_a[b↔ (b+1)]c_a(b+1) 
    	+	  k_a[(b-1)↔ b] c_a(b-1).

Here we have defined

    k_[(a-1)↔ a]b =
        k_0  exp(i∑_j=1^Λ        A^j_(a-1)b-A^j_ab/2 α_jsin (α_j z) )

and 

    k_a[(b-1)↔ b] =
        k_0  exp(i∑_j=1^Λ        A^j_a(b-1)-A^j_ab/2 α_jsin (α_j z) ).


Next we assume that the refractive index is oscillating fast. That is, the scale of variation in the z-direction of the modes c_ab(z) is much larger than 1/α_j for all j. 
Furthermore we choose

    α_j=q_j ω,

with q_j positive natural numbers. 
Integrating (<ref>) over one modulation period T:= 2π/ω under these assumptions then yields the effective coupling constant

    k_[(a-1)↔ a]b^eff :=1/T∫_z_0^z_0 +T k_[(a-1)↔ a]b(z)dz

and similarly in the b-direction. We now expand the integrand using the generating function of the Bessel functions of the first kind

    exp(i σsin(x) ) = ∑_n∈Z J_n(σ)e^inx,

for each index j. 
In Appendix <ref> we show that one may neglect contributions that contain J_n ≠ 0(...) terms.

Therefore, we may approximate our coupling constants as

    k_[(a-1)↔ a]b≈ k_0  ∏_j=1^Λ J_0(  A^j_(a-1)b-A^j_ab/2 α^j)

and

    k_a[(b-1)↔ b]≈ k_0  ∏_j=1^Λ J_0(  A^j_a(b-1)-A^j_ab/2 α^j).

Later, we will need to make the coupling between neighbouring waveguides position dependent. Suppose we are given a scalar function h on the lattice sites:

    0 < h(a,b) ≤max_c,d h(c,d)

For reasons that will become apparent soon, we would like to choose the coupling between sites to be proportional the geometric mean of the values given by the values the function h takes at the respective sites times the original coupling constant. This would for example require the coupling between c_ab and c_(a+1)b to equal k √(h(a,b)h((a+1),b))/max_c,d h(c,d).
To achieve this we introduce two refractive index modulations, whose amplitudes we will denote by A_ab and B_ab respectively and whose modulation speeds are denoted by α and β.
For the "A" modulations we assign amplitudes according to

    [
    0         A_(a+1)(b+1)
             
    
    A_ab        0   
    ],

extended periodically over the entire lattice and where A_ab is chosen such that

    J_0(A_ab/2α)=√(h(a,b)/max_c,d h(c,d)).

For the "B" modulation, we assign amplitudes according to 

    [
    B_a(b+1)        0  
             
    
    0         B_(a+1)b
    ],

with 

    J_0(B_ab/2β)=√(h(a,b)/max_c,d h(c,d)).

This then yields coupling constants according to 

    k_[a ↔ (a+1)]b =k   √(h(a,b)h((a+1),b))/max_c,d h(c,d)

and similarly in the b-direction.




§ WAVEGUIDE IMPLEMENTATION OF THE DIRAC EQUATION

In this section we describe how one may implement the Dirac equation in a coupled waveguide array. Let us first consider the case without gravity and electromagnetism.



 §.§ The flat, uncharged case


The Dirac equation for a free particle in 2+1 dimensional flat spacetime can be written as

    i ∂_t Ψ = -i (σ_y ∂_x + σ_x ∂_y )Ψ + m σ_z Ψ.

To bring this equation onto a square lattice, we will employ equations (<ref>) and (<ref>).
We denote the discretized two-spinor components as

    jψ_n^m:=Ψ_j(δ_x + nδ_x ,δ_y+mδ_y).

On a `unit cell', let us try assigning spinor components as 

    [ 1ψ_n-1^m+2    2ψ_n^m+2    1ψ_n^m+2    2ψ_n+1^m+2
    
    i 2ψ_n-1^m+1    i 1ψ_n^m+1    i 2ψ_n^m+1    i 1ψ_n+1^m+1
     1ψ_n-1^m    2ψ_n^m    1ψ_n^m    2ψ_n+1^m
    
    i 2ψ_n-1^m-1    i 1ψ_n^m-1    i 2ψ_n^m-1    i 1ψ_n+1^m-1
    ]

which is to be periodically extended over the entire lattice.
The derivatives are discretized as

    ∂_y Ψ_i ≈Ψ_i(n δ_x ,(m+1) δ_y) - Ψ_i(nδ_x ,(m-1)δ_y)/2δ_y

and

    ∂_x Ψ_i ≈Ψ_i((n+1) δ_x ,m δ_y) - Ψ_i(nδ_x ,mδ_y)/δ_x

or

    ∂_x Ψ_i ≈Ψ_i(n δ_x ,m δ_y) - Ψ_i((n-1)δ_x ,mδ_y)/δ_x,

depending on the spinor component "i" and the column of the lattice.

Now suppose that the coupling constants between any two given lattice sites always have the same absolute value. Then, to realize the discretized version of (<ref>), we necessarily have to set the absolute value of our coupling constant to

    |k| = 1/δ_x=1/2δ_y.

A quick calculation shows that the signs of the coupling constants have to vary as

    [ψ    -    ψ    +   ψ   -    ψ
    
       +         +            +           +         
    ψ    +    ψ    -   ψ   +    ψ
    
    -            -            -            -     
    ψ    -    ψ    +   ψ   -    ψ
    
    +            +            +           +      
    ψ    +    ψ    -   ψ   +    ψ
    ].

As it turns out, it is impossible to achieve such a prescription using only a single modulation of the refractive index. Since our effective coupling constants scales with the number of introduced modulations Λ as   [J_0(...)]^Λ, we want to minimize the number of introduced refractive index modulations.

To achieve this goal we note that multiplying an entry in (<ref>) with (-1) changes the sign of the immediately surrounding coupling constants that need to be achieved. This means that if we change (<ref>)  to

    [
    - 1ψ_n-1^m+2    2ψ_n^m+2    1ψ_n^m+2     - 2ψ_n+1^m+2
    
    
    i 2ψ_n-1^m+1    -i 1ψ_n^m+1    -i 2ψ_n^m+1    i 1ψ_n+1^m+1
     1ψ_n-1^m    - 2ψ_n^m    - 1ψ_n^m    2ψ_n+1^m
    
    -i 2ψ_n-1^m-1    i 1ψ_n^m-1    i 2ψ_n^m-1    -i 1ψ_n+1^m-1
    ],

(<ref>) changes to

    [ψ    +    ψ    +   ψ   +    ψ
    
    -         -           -          -      
    ψ    -    ψ    -   ψ   -    ψ
    
    -            -            -            -        
    ψ    +    ψ    +   ψ   +    ψ
    
    -           -            -           -       
    ψ    -    ψ    -   ψ   -    ψ
    ],
 
which can be achieved with only one sinusoidal modulation

    σ_mn(z)=C_mn/2cos(γ z).

To specify the required amplitude prescription, we note that the relation

    J_0(ξ) = -J_0(3ξ)

has a non-empty set of solutions. The first few positive roots are given by ξ_1 ≈  2.704, ξ_2 ≈  5.83 and ξ_3 ≈ 8.97 and e.g. J_0(ξ_1)=-J_0(3ξ_1 ) ≈ -0.144.
We define  C:= ξ_1/2γ and assign the amplitudes as

    [
    	0     +     3C     +   0    +    3C   
    
    	-         -           -          -        
    
    	C     -     2C     -   C    -    2C     
    
    	-            -            -        -         
    
    		0     +     3C     +   0    +    3C   
    
    	-           -            -        -        
    
    	C     -     2C     -   C    -    2C   
    ].
 
where the signs of the coupling constants follow from the assigned amplitudes by equations (<ref>) and (<ref>).

We have thus realized the discretized massless 2 + 1 dimensional Dirac equation with

    J_0(3ξ_1)=|k| = 1/δ_x=1/2δ_y.

If we would like to take δ_x =  δ_y, so that the discretisation lengths in the two directions coincide, it is obvious from the above considerations that this corresponds to introducing different coupling constants in the x and y directions:

    k_x=1/2 k_y.

This can be achieved by changing the distance of the waveguides in an experimental setup, or by introducing additional refractive index modulation in the x direction similarly for every y column, according to [0,E] extended periodically so that J_0(E/2 ϵ) = 1/2. This changes the effective coupling in the x direction to one half of its original size, while leaving the coupling in the y direction invariant.
We will however not do so and instead work with different discretisation lengths in the x and y directions, so that the coupling constants in both directions can have the same value.
Including a mass term can be achieved through the usual method of alternating high and shallow refractive indices <cit.>.




 §.§ The general case



Our starting point is (<ref>), which we repeat here for convenience:

    ∂_tχ=   -ie^ρ(σ_y∂_x +σ_x∂_y )χ
       +(σ_ze^Φ m + 1[ϕ- ∂_t Λ] ) χ

Just like the mass term, (σ_ze^Φ m + 1[ϕ- ∂_t Λ] ) can be easily implemented by changing the refractive index of the respective waveguide appropriately.
Hence we turn to the derivative term.
For ease in notation we define f:= e^ρ. 
The relevant part of the equation may thus be written as

    i∂_tχ= -if(σ_y  ∂_x + σ_x ∂_y  ) χ

If we want to put this on the lattice, we will have to make the absolute values of the effective coupling constants position dependent. Furthermore, we need to be careful to obtain the correct differential equation in the limit of vanishing discretization length.
For this, the effective couplings of the spinor components (as opposed to the amplitudes in waveguides) to its neighbours are the determining factors.

It is easy to verify that in the discretized version of (<ref>) the coupling of a spinor component ψ_n^m to its neighbours jψ_k^l must be local, i.e. the absolute value may only depend on the tuple (nm) and must not depend on the coordinates of the neighbouring sites (kl). 

To achieve this, we will use the lattice constructed in the previous subsection, superimposed with the modulation of coupling constants that we introduced at the end of section <ref>.
Note that while we were considering a function h(a,b) on the lattice sites {(ab)} there, we are now considering a function f(n,m) on discretized space points {(nm)}. Given the prescription (<ref>) of spinor components, the two are related as

    h(a,b) := f(⌊  a/2 ⌋,b),

as is easily seen.
Here ⌊·⌋ denotes the floor function that assigns to each real number its closest smaller integer.


Additionally we modify the assignment of spinor components to our lattice by  pointwise multiplication with 

    [
    Z_(n-1)(m+2)    Z_n(m+2)    Z_n(m+2)     Z_(n+1)(m+2)   
    
    Z_(n-1)(m+1)    Z_n(m+1)    Z_n(m+1)    Z_(n+1)(m+1)   
    
    Z_(n-1)m    Z_nm    Z_nm    Z_(n+1)m   
    
    Z_(n-1)(m-1)    Z_n(m-1)     Z_n(m-1)    Z_(n+1)(m-1)   
    ]

extended periodically.
We now define

    √(f(n,m)/f(n+1,m)) Z_nm = Z_(n+1)m

and likewise

    √(f(nm)/f(n,m+1)) Z_nm = Z_n(m+1).

The consistency of these conditions is easily verified.
We have one degree of freedom left by choosing the initial Z value in the inductive definition above.
Let us fix 

    Z_00=1/√(f(0,0)).

Then we have

    Z_mn   =1/√(f(0,0))·∏_i=1^n  √(f(i-1,0)/f(i,0))·∏_j=1^n √(f(n,j-1)/f(n,j))
       =1/√(f(n,m)).

Note that we could have chosen any other path from (0,0) to (nm).
Plugging in the relations  (<ref>) -  (<ref>) into the effective coupled mode equations that we constructed in Section <ref> then shows that we have indeed realized the discretized version of equation (<ref>) if we choose the initial coupling constant k_0 according to 

    k_0:=  max_m,n f(n,m)/δ_x = max_m,n f(n,m)/2δ_y.


Checking that the discretisation scheme yields the correct differential equation in the limit of vanishing stepsize is straightforward, since the function multiplying the discretized derivative term is local.
Lastly we note that our implementation scheme is able to precisely realize equations of type (<ref>). The only part that is not immediately obvious is 
the relationship between the functions multiplying the spatial derivative terms. 
To unravel this relationship, assume that the x-derivative is multiplied by f=e^q(x,y) and the y derivative is multiplied by g=e^p(x,y).
For consistency we must then require

    g_n,m/g_n,m+1f_n,m+1/f_n+1,m+1=f_n,m/f_n+1,mg_n+1,m/g_n+1,m+1

for arbitrary δ_x,δ_y.
Taking the logarithm of this equation, dividing by (δ_x  δ_y) and taking the limits δ_x,δ_y → 0 then implies

    ∂_x ∂_y q = ∂_x ∂_y p

by Taylor's theorem.
Hence

    p - q =  B(x) + C(y) ,

for some functions B and C of only a single variable.
Thus the partial differential equations that are implementable indeed take the form of (<ref>).


Now that we have shown, in principle, how the Dirac equation in classical backgrounds may be implemented in waveguide arrays, we would like to draw the reader's attention to the considerable engineering challenge associated with actually fabricating the waveguide arrays needed for our implementation scheme:
In addition to the non trivial coupling between the spinor components, we have encoded the entire information about the mass of our test particle, the background spacetime and the electromagnetic background it experiences into the refractive indices of the waveguides. The limiting case of a flat spacetime and the gauge field set to zero is certainly realisable, as there is only a single sinusoidal refractive index modulation required (for an experimental realisation of a similar setup see e.g. <cit.>).
It has to be investigated on a case by case basis for any given more complex background scenario if it might be realised with state of the art techniques.



§ GRAVITATIONAL AHARONOV-BOHM EFFECT


In this section, we devise a thought experiment to observe the gravitational Aharonov-Bohm effect. We then show that an analogue can be implemented in a tabletop experiment with the methods constructed above. 


 §.§ Theoretical considerations

Our starting point is the conical space. As argued in Appendix <ref>, we may write the metric as

    ds^2    = dt^2 - (√(x^2+y^2))^-2 Δ(dx^2+dy^2),

in a global coordinate system.
Going to polar coordinates, we can write this as

    ds^2 = dt^2-R^-2Δ(dR^2+R^2 d τ^2).

The spacetime is locally flat, i.e. every point has an open coordinate neighbourhood on which the metric takes Minkowskian form. While there is no such coordinate chart that covers all of the cone, one may however introduce coordinate neighbourhoods in which the metric is Minkowskian, that cover the entire cone except for a straight line starting at the tip of the cone. For fixed time, this can be understood geometrically as cutting the cone along this line and flattening out the obtained surface <cit.>. 

The result is a copy of R^2 with a wedge of angular size Δ2π removed. In such a coordinate system, the metric is given by 

    ds^2 = dt^2-dρ^2 - ρ^2 dθ^2

with

    0 < θ < (1-Δ) 2 π.

The points for which θ = 0 and θ = (1-Δ) 2 π, i.e. at the edge of the removed wedge, would correspond to the line along which we cut the cone and need to be excluded. This is for the same reason one excludes the values {0, 2π} in the usual polar coordinates of R^2. To be mathematically precise, this ensures that our coordinate maps are  diffeomorphisms between open sets. Thus we see that to describe the entire cone using only flat coordinate systems, we need an additional flat coordinate patch. This corresponds to cutting the cone along a different line. To signify the difference between these coordinate systems, we denote the angular variable by ϕ for the second coordinate system. As before, we have

    ds^2 = dt^2-dρ^2 - ρ^2 dϕ^2

with

    0 < ϕ < (1-Δ) 2 π.


We make the cuts along the positive and negative x axes of our global coordinate system (see Figure <ref>). 

The idea now is to mimic the Aharonov-Bohm effect for U(1) gauge fields. For this we want to create an input state that corresponds to two  separately localised parts Ψ_0, Φ_0  of a wavefunction with initial momenta chosen such that one of the localised parts will travel above the tip of the cone at (x=0,y=0), while the other travels below. Heuristically, we think of the two areas of localisation of the wavefunction as two separate particles

In flat spacetime it is well known how to create  localised  states with  Gaussian momentum distribution. One e.g. projects a plane wave solution in momentum space onto the positive energy subspace, multiplies with a Gaussian centered around the desired momentum and Fourier transforms <cit.>.
We understand the free time evolution of such initial states in the flat spacetime and can easily calculate it numerically. We can use this knowledge to understand the dynamics in conical space, by making use of the fact that our spacetime is locally flat.
We initially work in the coordinate system characterised by the angular variable θ and prepare an input state

    Ψ_flat := Ψ_0 + Φ_0

normalised as ||Ψ_0||=||Φ_0||=1/√(2).
 Ψ_0 is taken to be localised in sector II and Φ_0 in sector III (compare figure <ref>).

We choose the θ-coordinate system for the preparation of the input states, since we heuristically think about the two localised parts of the initial state as originating from a common source located in the direction opposite to the initial momenta of the two parts. Hence we need a coordinate patch that covers the entire geometry of such a setup, including the point of the source, the parts of the spacetime where the states Ψ_0,Φ_0 are localised and the trajectory between the source and the state. We want to describe the subsequent time evolution in the coordinate system characterised by the angular variable ϕ, since we will mainly be concerned with observations in the sectors I and IV and hence want a flat coordinate system that is able to describe this patch. 


To be able to do so, we require the time evolved states to not extend through the line along which we cut our cone to go to the ϕ-coordinate system (i.e. {x<0, y=0} under the identification of the constant time slice of the manifold with the global coordinate system). 
Fortunately enough this can be guaranteed:
Let us first approximate the initial Gaussian wavefunction by zero whenever the values fall below a certain threshold. This can be done smoothly, if we subsequently mollify around the cutoff. Since the time evolution in flat space is implemented as a unitary operator, the time evolved truncated state approximates the time evolved true state with the initial error in norm at all times on the flat patch. Furthermore we note that the Dirac equation imposes a finite propagation speed. That is, for continuously differentiable ψ∈ C^1(R^2) which is non-zero only within a ball of radius r centered at x_0 (i.e. suppψ(0,·) ⊂ B_r(x_0)), we have suppψ(t,·) ⊂ B_r+t (x_0), i.e. it is non-zero only within a Ball of radius (r+t) <cit.>. 
Hence if we prepare the initial states around large enough initial momenta and far enough away from the half line {x<0, y=0}, we avoid the truncated states to reach this half line  before the expectation values have propagated beyond the tip of the cone in this coordinate system.

Hence we may change our coordinate system to the one characterized by the angular variable ϕ and do the subsequent calculations in this coordinate system.
We have to first check how the input states change under such a coordinate transformation to find their form in the new coordinates. 
For this, we need to know the angular change ϕ - θ.
We obtain this information by making reference to the angular variable τ of the global coordinate system.
From the relation of metrics in the respective coordinate systems (cf. Appendix <ref>), we can infer that where the coordinate systems overlap

    dϕ = (1-Δ) d τ = d θ .

Furthermore the half line

    {lim_τ↓ 0^+R(cos(τ),sin(τ))^T | R>0 }

in the coordinate system characterized by τ corresponds to   

    {lim_θ↓ 0^+ρ(cos(θ)sin(θ))^T  | ρ > 0 }

in the flat coordinate system characterized by θ (see Figure <ref>). Here

    ρ = (1-Δ)^-1R^1-Δ.


Strictly speaking we need to take the limits on the manifold itself and not in the coordinate representation, since as argued above, the limits of the considered sequences are not in the coordinate images anymore. We however stick to this notation since it is less cluttered and there is no room for confusion.

For y<0 we know that the half line 

    {R(cos(π),sin(π))^T | R>0 }

in the coordinate system characterized by τ corresponds to the half line

    {lim_ϕ↓ 0^+ρ(cos(ϕ),sin(ϕ))^T  | ρ > 0 }.

Integrating the differential equations (<ref>) with these initial conditions then yields

    θ   = (1-Δ)τ, 
    ϕ   =(1-Δ)(τ-π).

For the set where y>0, we know that 

    {R(cos(π),sin(π))^T | R>0 }

corresponds to 

    {lim_ϕ↑ (1-Δ)2πρ(cos(ϕ)sin(ϕ))^T  | ρ > 0 }.

Together with (<ref>), this then implies

    θ   =(1-Δ)τ
    ϕ   =(1 - Δ)(τ +  π).


An easy calculation then shows

    [ϕ - θ]_mod2π = (1-Δ)π =: α

for y> 0, and

    [ϕ - θ]_mod2π =(1+Δ)π =: β

for y< 0.
Of course we can infer the value of the change in the angular variable only up to 2 π. 
Since we are working with spinors, this is not a triviality. As we will see momentarily, it does not pose an obstacle though.

Under a Lorentz transformation Λ, infinitesimally characterized by ω_ab(Λ), spinors change according to <cit.>

    ψ(x) ⟼Λ_1/2(ω_ab(Λ))ψ(Λ^-1x).

Here, we are working with the spin-1/2 representation, denoted by Λ_1/2(·).
Given the Clifford algebra structure of {γ_0,γ_1,γ_2}, we can represent the effect of a Lorentz transformation Λ through <cit.>

    Λ_1/2(ω_ab(Λ) ) := exp( i/2ω_ab(Λ)i/4[γ^a,γ^b]).

In our case of a rotation Λ = R_2(η) in the plane, the only non-vanishing entries of ω_ab are ω_12=-ω_21= η. This then yields 

    V(η) :=Λ_1/2(ω_ab(η)) = exp( i/2ησ_z).

The two coordinate systems dissect the underlying spatial slice of our manifold into two open disjoint sub-spaces, the one where y >0 and the one where y<0. Thus we a priori do not know to which of the two possible elements {V(η),V(η + 2 π)) } the elements  R_2(η) ∈ SO(2), where η = α,  β, need to be lifted for a consistent prescription. We can however infer this by considering the case Δ→ 0 for which the transformation for η=α needs to coincide with the transformation for η = β. As this is the case for the choice of representations of the angular variables  in (<ref>) and (<ref>) we see that we are consistent in simply inserting η=α, β into (<ref>).
We therefore infer that our input states must change according to

    Ψ_0(x) ⟼  V((1-Δ)π)Ψ_0(R_2^-1((1-Δ)π)x)
       Φ_0(x) ⟼ V((1+Δ)π) Φ_0(R_2^-1((1+Δ)π)x).


Next, let us see how to describe the time evolution of the input states that we have in the coordinate system with angular variable ϕ.
As was to be expected, we find that it does not matter if one first time evolves and then changes the coordinate system, or if one first changes the coordinate system and then time evolves in the new coordinate system.
To see this rigorously, let us consider a Lorentz transformation Λ: R^1,2→R^1,2, mapping from the coordinates {y'^μ} to {z^μ}. Recalling the covariance of the flat Dirac equation, we know that

    (i γ^μ(Λ^-1)^ν_μ∂/∂ z^ν-m)Λ_1/2(ω_ab(Λ)) ψ(t,Λ^-1z)
    
    = Λ_1/2(ω_ab(Λ)) (i γ^μ∂/∂ y'^μ-m) ψ(t,y') =0

holds for solutions ψ(t,x) <cit.>. In our case, t=z^0=y'^0 and the transformation Λ is a pure rotation R_2(η):R^2 →R^2.
Since in this case [ γ_0 , Λ_1/2(ω_ab(η))]=0, which in our choice of representation is easily seen since γ_0 = σ_z and Λ_1/2(ω_ab(R_2(η))) = e^i/2ησ_z, we have for the time evolution U_z(t) in the z-coordinate system

    U_z(t)V(η) Ψ_0(R_2^-1(η)z)
    = V(η) Ψ(t,R_2(-η)z).

where Ψ(t,y') = U_y'(t) Ψ_0(y') is the solution of the Dirac equation in the y'-coordinate system. 

Now first take two copies of R^2 with coordinates

    (y'^1,y'^2)    := r(cos(θ),sin(θ))  : 0< θ < 2 π
    
    (z^1,z^2)    := r(cos(ϕ),sin(ϕ))  :  0< ϕ < 2 π,

Then we can think of the coordinate representations that describe parts of our manifold (i.e. the ones where 0< ϕ,θ < (1-Δ) 2 π)) as the open subsets of these copies of R^2 in which the corresponding angular wedge is removed.
If we take η = α, we can completely describe the time evolution of the initial Ψ_0 in this setting:

    U_z(t)V(α) Ψ_0(R_2^-1(α)z)
    = V(α) Ψ(t,R_2(-α)z).

We only have to make sure that at no point, after transforming the state to the z-coordinate system are there non-zero contributions in the wedge (1-Δ)2 π≤ϕ≤ 2 π in order to maintain a consistent description.
We argue similarly for the case of the initial Φ_0; only here we have to set η=β:

    U_z(t)V(β) Φ_0(R^-1(β)z)
    = V(β) Φ(t,R_2(-β)z).


In total, assuming that neither state is non-zero in the missing wedge in the z-coordinate system, we have as our total state at time t

    Ψ_tot(t,z) = V(α) Ψ(t,R_2(-α)z) +V(β) Φ(t,R_2(-β)z).

The probability density everywhere except on the half line { y=0,x<0 } and thus especially in the quadrants I and IV, expressed in coordinates z^1=rcos(ϕ), z^2=rsin(ϕ) is then given by

    ρ_Ψ+Φ(t,z)= Ψ̅Ψ(t,R_2(-α)z)+Φ̅Φ(t,R_2(-β)z)
        + 2{Ψ̅(t,R_2(-α)z) e^iσ_zΔπΦ(t,R_2(-β)z)  }.

This dependence on the deficit angle is precisely what we expected after our discussion in section <ref>. 

We note the similarity to the corresponding result for the interference pattern of electrons under the influence of the electromagnetic Aharonov-Bohm effect. There, the role of the deficit angle is played by the magnetic flux:
For the Aharonov-Bohm effect, one can obtain time-dependent solutions by time evolving initial states as if no gauge potential is present as long as we don't have to consider paths encircling the solenoid containing the magnetic field (i.e. in any region of the space with vanishing fundamental group). Then, to describe the interference pattern of parts of the wavefunction that have `travelled' above and below the solenoid, we need to gauge transform to a common gauge choice. Since the gauge fields can't be set to zero globally, this  modifies the time evolved states by adding a phase corresponding to whether the states travelled above or below the solenoid <cit.>.

In the gravitational case, we time evolve initial states as if our flat coordinate system extends indefinitely. This description is valid as long as we don't have to consider paths that encircle the tip of the cone. 
To calculate the probability density on the side of the cone opposite to where we prepared the initial state, we go to the coordinate system that completely covers this side of the cone.

Spinors transform non-trivially under rotations; not only 
is their argument transformed as x ↦ R^-1(η)x but they are also multiplied with an element of Spin(2). In our representation, this is given by 

    V(η) = e^iη/2σ_z

 and can thus heuristically be understood as a generalized phase factor.
 Depending on whether they travelled above or below the cone, the respective localized parts of the wavefunction are multiplied with different generalized phase factors, according to η = α or η = β.

We note that with our coupled waveguide array tuned to host a conical geometry, we can realize such an interference experiment directly, since the intensity in the lattice can be used to infer the probability density, which transforms as a scalar under a general coordinate change <cit.>. 
We now want to detail how to achieve this.



 §.§ Physical implementation in a waveguide array

Our starting point is (<ref>), repeated here for convenience:

    i ∂_t χ =i f( σ_y ∂_x + σ_x ∂_y) χ+ mσ_z χ

Here 

    f = 1/g = (√(x^2+y^2))^Δ,

and χ is a rescaled spinor, related to the true spinor as 

    χ:= √(g)ψ = 1/√(f)ψ.

We discretize the space as

    [ x; y ]⟼[   (1/2+n)δ_x; ( 1/2 +m)δ_y ]  m,n ∈Z,

which is useful because

    lim_x,y → 0 g(x,y) ⟶∞,

and we can cap this potentially singular behaviour in the field redefinition (<ref>) using (<ref>) since no space point is directly placed in the center.
We define

    jχ_n^m:=χ_j(( 1/2 + n) δ_x ,( 1/2 + m)δ_y),

as well as

    f(n,m):=   (( ( 1/2 + n) δ_x)^2+( ( 1/2 + m) δ_y)^2)^Δ/2,


    √(g(n,m)) := 1/f(n,m)

and

    h(a,b) := f(⌊ a/2 ⌋,b).


We then take our waveguide lattice with three distinct modulations of the refractive index.
We assign rescaled spinor components jχ_n^m according to (<ref>), further rescaled by (<ref>). 
To fix the values {Z_mn}, let us choose

    Z_00=1/√(f(0,0)).

Then we have by (<ref>)

    Z_mn=1/√(f(m,n)) =√(g(n,m)).

Denoting

    (g jψ)_n^m := g(n,m) · jψ_n^m,

we are thus assigning rescaled spinor components according to:

    [
    -(g 1ψ)_n-1^m+2   (g 2ψ)_n^m+2   (g 1ψ)_n^m+2     -(g 2ψ)_n+1^m+2
    
    
    i(g 2ψ)_n-1^m+1    -i(g 1ψ)_n^m+1    -i(g 2ψ)_n^m+1    i(g 1ψ)_n+1^m+1
    (g 1ψ)_n-1^m    -(g 2ψ)_n^m    -(g 1ψ)_n^m   (g 2ψ)_n+1^m
    
    
    -i(g 2ψ)_n-1^m-1    i(g 1ψ)_n^m-1    i(g 2ψ)_n^m-1    -i(g 1ψ)_n+1^m-1
    ]



To achieve the correct signs of the effective coupling constants we employ an amplitude prescription according to (<ref>) for the first modulation. 
The two remaining modulations are then used to achieve the necessary position dependent modulation of the absolute value of the coupling constants.
As detailed in (<ref>) and (<ref>), every site that has a non-zero amplitude for  one modulation has zero amplitude for the other and vice versa. Furthermore every non-zero amplitude of a given modulation is surrounded by sites that have zero amplitude for this modulation.
In the present case, the amplitudes on the lattice site (ab) are then chosen such that 

    J_0(Ξ_ab) = √(h(a,b)/max_n,mf(n,m)),

where in the notation of section <ref>, Ξ_ab=A_ab/α or Ξ_ab=B_ab/β depending on which amplitude is non-zero and we the maximum  is taken over all {(mn)} considered in the actual setup.
This then yields the correct coupling constants for each wave guide (a,b) to its neighbours to ensure the correct coupling of the spinor component jΨ_ n^m to its neighbours.

Next let us consider the input state preparation.
Given localized initial states Ψ_0(y') and  Φ_0(y') in the coordinate system characterized by the angular variable θ, we need to know how to express them in the global coordinate system characterized by τ.
We know

    (1-Δ)τ = θ,

and hence

    [τ - θ]_mod2π = Δτ.

Since the vielbein is diagonal, the combined initial state Ψ_flat(y') in the flat θ-coordinate system then simply corresponds to the state

    Ψ_ initial= V(Δτ(x,y))Ψ_flat(y'(x,y)) .

Here τ≡τ(x,y) is given by the relation

    (x,y)= √(x^2+y^2)(cos(τ(x,y)),sin(τ(x,y))),

and (cf. (<ref>))

    (y'^1,y'^2)(x,y)= √((y'^1)^2 + (y'^2)^2) (cos(θ), sin(θ))
       =√(x^2+y^2)^1-Δ/(1-Δ) (cos((1-Δ)τ), sin((1-Δ)τ)).


As input state in our lattice, we then prepare the rescaled discretized verion of of (<ref>), given by

    j(Ψ_disc.)_ n^m:= g(n.m)    j(Ψ_initial)_ n^m.

The evolution along the waveguides then mimics the time evolution of the approximated quantum system, i.e. that of a "particle" governed by the Dirac equation in conical space, in the global coordinate system. Hence if we want to know the probability density of the quantum  system at a given time, we can infer it from the intensity 

    I_ab= |c_ab|^2

 in our waveguide system at the corresponding transverse distance from the points where the input state was prepared.

To observe the "phase shift" of the gravitational Aharonov-Bohm effect, let us now prepare the input state

    Ψ_flat(y') := Ψ_0(y')+ Φ_0(y'),

as discussed above in section <ref> (c.f. also Figure <ref>).
Note that the demands put on the  waveguide system that needs to be fabricated by the requirements on the position of the localised parts of the initial state can be mitigated by varying the discretization length or the physical (non-effective) coupling constants.  The requirement of the input state  being predominantly composed of modes corresponding to large enough momenta however, leads to a potentially complicated phase structure that might vary rapidly with the waveguide position for the input state. This might add additional challenges to the optical input state preparation.

After the input state is prepared and some "time" has passed, the probability density will then be given by (<ref>) with z=z(x,y), since it simply transforms as a scalar under coordinate changes <cit.>.
Note that 

    (z^1,z^2)(x,y) = r(cos(ϕ),sin(ϕ))
       =√(x^2+y^2)^1-Δ/(1-Δ)(cos(ϕ(τ(x,y))),sin(ϕ(τ(x,y)))).

The relation ϕ(τ) between the angular coordinates is given by (<ref>) and (<ref>).
These equations imply 

    lim_τ↓ 0ϕ(τ) = (1-Δ)π = lim_τ↑ 2πϕ(τ),

and thus ϕ(τ) is well defined.

In our waveguide array experiment, the amplitudes in the respective waveguides correspond to the spatially discretized version of the rescaled spinor components. Hence we can reconstruct the discretized probability density as 

    ρ_Ψ+Φ(t,z([1/2 + n ]δ_x,[1/2 + m ]δ_y))
       =1/(g(n,m))^2[I_(2n)m + I_(2n+1)m]
       = (f(n,m))^2 [I_(2n)m + I_(2n+1)m].

Obtaining these experimentally and subsequently plotting them in the z coordinate system (i.e. as a function of z and not as a function of (x,y)) using the relation z(x,y), then precisely yields a discretized approximation of (<ref>).

Hence we can observe the deficit angle as the Δ-dependent shift of this interference pattern. One might obtain the numerical value of Δ from the interference experiment described above by numerically calculating or obtaining from a separate waveguide experiment tuned to host flat space, the time evolution of the initial states Ψ_0, Φ_0 in the flat coordinate system. Then the function

    f_Δ(z):=  Ψ̅Ψ(t,R_2(-(1-Δ))z)
       	+Φ̅Φ(t,R_2(-(1+Δ)z)
        + 2{Ψ̅(t,R_2(-(1-Δ))z) e^iσ_zΔπΦ(t,R_2(-(1+Δ))z)  }.

at fixed time t can be fitted to the data obtained in the waveguide experiment detailed above. Here the fit proceeds with parameter Δ. The resulting value of Δ yields the value of the "phase shift" e^iΔπ. 




§ CONCLUSION

We developed a scheme to implement the Dirac equation in 2+1 dimensions in curved spacetime and classical electromagnetic background in a coupled waveguide array.
To do so, we had to introduce effective, site dependent coupling constants. This was achieved by periodically varying the refractive index of the respective waveguides. In the case of flat spacetime, one modulation sufficed, since only varying signs of the coupling constants needed to be achieved. For the general case, two more modulations and a spinor rescaling enabled the introduction of arbitrary site dependent couplings of spinor components. The electromagnetic potential was encoded into the refractive index too.
The implementation scheme in principle allows for implementations of many physically interesting settings. As an example, we devised a thought interference-experiment to observe a gravitational analogue of the Aharonov-Bohm effect in the spacetime of a cone. The global geometry of this locally flat setting introduced a generalized phase difference  between the time evolution of separately localized states.
Finally, we briefly explained how an analogue of this setting might in principle be implemented and how the associated generalized phase shift could be observed in a waveguide array.

 

Acknowledgements: 
C.K. acknowledges travel support provided by the Elite Master Program Theoretical and Mathematical Physics of the Free State of Bavaria, Germany. 
	D.A. was supported by the National Research Foundation, Prime
Minister’s Office, Singapore and the Ministry of Education, Singapore
under the Research Centres of Excellence programme. D.A. was also
partially funded by the Polisimulator project co-financed by Greece and
the EU Regional Development Fund, the European Research Council under
the European Union’s Seventh Framework Programme (FP7/2007–2013).






tocsectionAppendices





§ APPENDICES




 §.§ Introduction to Gravity in 2+1 dimensions

Since certain aspects are quite different from the more familiar 3+1 dimensional case, it is expedient to briefly discuss the most distinctive features of three dimensional gravity. We closely follow <cit.>.


  §.§.§ General remarks on 2+1 dimensional gravity

Our starting point is a three dimensional pseudo Riemannian manifold (M,g). One can show that in 2+1 dimensions, the Ricci tensor completely determines the Riemann tensor <cit.>:

    R_μνρσ   = g_μρ R_νσ + g_νσ R_μρ - g_νρ R_μσ - g_μσ R_νρ
       -1/2(g_μρg_νσ-g_μσg_νρ)R.

Gravity is defined via Einstein's field equations

    G_μν=8π T_μν,

where we used the definition of the Einstein tensor,

    G_μν := R_μν - 1/2 R g_μν,

and absorbed a possible cosmological constant into the Stress-energy tensor.
Note that throughout this paper we work in Planck units.
Since the vanishing of Einstein's Tensor implies a vanishing Ricci salar, it also implies that the Ricci tensor vanishes.
But in 2+1 dimensions, this now immediately gives a vanishing Riemann tensor. 
There is thus no way that matter can communicate gravitationally in three dimensional spacetime and curvature at a specific point of the manifold is equivalent to matter being present at this point.
All effects that localized sources have are on the global geometry <cit.>. 

Another interesting feature for this dimensionality is the fact that the Cotton tensor, given in coordinates by

    C_ijk:=∇_k R_ij-∇jR_ik+ 1/4( ∇_j Rg_ik   -∇_kRg_ij),

plays the role of the Weyl tensor of higher dimensional pseudo Riemannian manifolds: In 2+1 dimensions, the vanishing of the Cotton tensor is equivalent to conformal flatness, while in higher dimensions, conformal flatness is equivalent to the vanishing of the Weyl tensor.
It is not hard to see that a vanishing Ricci tensor implies a vanishing Cotton tensor <cit.>. Thus vacuum solutions to Einstein's field equations yield conformally flat metrics.
We now turn to a specific spacetime we want to study in greater detail.



  §.§.§ Example: The Cosmic string solution: Conical space

In 3+1 dimensions, the simplest cosmic string solution is a solution that is invariant under translations in one spatial direction <cit.>. It is thus essentially a 2+1 dimensional system with a trivial dependence on a fourth spacetime coordinate, to be called z in what follows. For reasons of generality we first discuss it in the context of 3+1 dimensional gravity. Then we will discard the z-coordinate and consider the 2+1 dimensional system. 

During a symmetry breaking phase, for example in the very early universe, topological defects may arise naturally <cit.>.
Cosmic strings are such hypothetical, almost-one-dimensional topological defects in spacetime, which are analoguous to other linear defects that are familiar from condensed matter systems <cit.>.
The system of equations governing the dynamics of the gravitating string can be made tractable by assuming a vanishing string diameter, thus approximating it as a line of zero width through a distributional δ-function energy-momentum tensor. Furthermore, one may assume a sufficiently weak gravitational field, so that the linearized Einstein equations provide a good approximation. For a detailed derivation, the reader is referred to <cit.> or <cit.> which we are following closely, below.

We consider a weak gravitational field in which the spacetime is almost Minkowskian,

    g_μν = η_μν + h_μν,

with |h_μν| ≪ 1. It is thus justified to linearize Einstein's equations in h_μν. Using the  harmonic gauge, specified by

    ∂_ν(h^ν_μ - 1/2δ^ν_μ  h^σ_σ) = 0,

one may find the solution of Einstein's equations for a straight string along the z axis, for which

    T^ν_ρ = μδ(x)δ(y)diag(1,0,0,1) .

The solution is found to be

    h_00=h_33=0 
       h:= h_11=h_22=μln(r/r_0)

with r:=√(x^2+y^2) and r_0 a constant of integration.

Consistency with the condition  |h_μν| ≪ 1 then limits the domain of applicability of the weak field equations. However, the current coordinate system is not a suitable one for this purpose. Instead we go to the cylindrical coordinates, in which

    ds^2=dt^2-dz^2-(1-h)(dr^2+r^2dθ^2).

Introducing 

    (1-8μ)ρ^2=(1-h)r^2

one obtains to linear order in μ,

    ds^2 = dt^2-dz^2-dρ^2 - (1-8μ)ρ^2 dθ^2.

Introducing a new angular coordinate 

    θ̃ = (1-4μ)θ

the metric can be written in the Minkowskian form

    ds^2 = dt^2-dz^2-dρ^2 -ρ^2 d θ̃^2

with

    0 < θ̃ < 2π (1-4μ).

again to linear order in μ.

One therefore sees that to first order in μ, the presence of the string introduces an azimuthal deficit angle 2 πΔ, where Δ≡ 4  μ.
while the spacetime remains locally flat.
The surface of constant t and z thus has the geometry of a cone.
The weak field approximation (with |h_μν| ≪ 1 ) is justified for μ≪ 1. Observations constrain G μ≲ 10^-5 for cosmological scenarios. 
The weak field metric   (<ref>) coincides with

    ds^2 = dt^2-dz^2-dρ^2 - (1-Δ)^2ρ^2 dθ^2

to first order.
The validity of this metric can be extended beyond linear perturbation theory by taking the spacetime metric to be given by (<ref>) outside the domain of applicability of linear perturbation theory <cit.>. 

The metric (<ref>) is locally Minkowskian and is thus an exact solution of the vacuum Einstein equations.  However it induces  a non-trivial global geometry. 
We will now focus on the study of the 2+1 dimensional version of this metric, which is obtained by discarding the dz^2-term. 
At this point one might wonder which 2+1 dimensional object could play the role of a source for this metric.
It turns out that the metric for a spinning point source in 2+1 dimensional gravity and conformal coordinates, is given by <cit.>

    ds^2= (dt + 4Jdθ)^2 - R^-8μ(dR^2+R^2 dθ^2),

with μ,J = const. and √(g)T^00=μδ(x_1)δ(x_2), T^ij=0, √(g)T^0i= 1/2Jϵ^ij∂_j δ(x_1)δ(x_2).
Here, x_1,x_2 are the cartesian coordinates corresponding to the cylindrical coordinate system the metric is expressed in.
Introducing a new radial coordinate 

    r:= (1-4μ)^-1R^(1-4μ),

this can be written as

    ds^2= (dt + 4Jdθ)^2 - dr^2 - (1-4μ)^2r^2 dθ^2.

It is then easy to see that the metric is flat outside the source.  As we have discussed before, this has to be the case in 2+1 dimensions.  
The spacetime outside an arbitrary distribution of particles in the centre of mass frame confined to a region r  r_0 can be described by this metric for r  r_0  <cit.>. Within this region the curvature is not necessarily zero as there are sources present.

It is not hard to see that for  J  0 the metric admits closed timelike curves (for example take t:= t_0, r:= r_0  4GJ(1-4gμ) and vary θ from 0 to 2π). To avoid this complication, we set J:= 0. Upon setting Δ:= 4μ we are thus left with metric (<ref>) after discarding the z-direction:

    ds^2 = dt^2-dρ^2 - (1-Δ)^2ρ^2 dθ^2.

This metric describes a conical space, which is a locally flat space with a wedge of angular size 2πΔ removed and the two lines, along which the wedge was cut out, identified.
Note that local flatness means that each point has an open coordinate environment on which the metric in coordinates is the Minkowski metric. This is different from simply introducing normal coordinates, in which case the metric a priori only takes Minkowski form at a single point.
The metric in a global coordinate system can be inferred from (<ref>) with J=0

    ds^2    = dt^2 - R^-2 Δ(dR^2+R^2 dθ^2)
       =dt^2 - (√(x^2+y^2))^-2 Δ(dx^2+dy^2).
 
This is clearly of the form (<ref>) with 

    e^Ψ = (√(x^2+y^2))^- Δ.

Since the spacetime is locally flat, it is clear that a test particle initially at rest relative to the string will stay at rest and there is no Newtonian gravitational force.
The global geometrical difference with the Minkowski space however gives rise to interesting non-trivial effects, like for example double images <cit.>. 

The situation is reminiscent of the Aharonov-Bohm effect in electrodynamics. There, one considers a punctured plane and a vector potential A that may be gauged away in any region with vanishing fundamental group, but not globally.
For a string, the curvature is confined to its core. We have chosen a vanishing diameter for simplicity, for the general case see <cit.>. The non-trivial global geometry however influences particles propagating in the flat region outside the core. A Minkowskian coordinate system can be chosen locally, but not globally.



 §.§ The Dirac equation in curved spacetime

In local coordinates, the Dirac equation is given by <cit.>

    [ iγ^μ∇_μ -m]ψ(x) = 0.

Here we are using the vielbein e^μ_a(x), defined by

    e^μ_a e^ν_bη^ab = g^μν,

with the Minkowski metric η to transform the local gamma matrices γ^a according to

    γ^μ(x) = e^μ_a(x)γ^a.

Note that {γ^μ(x),γ^ν(x) } = 2g^μν(x). 
The covariant derivative is given by

    ∇_ν :=  ∂_ν + Ω_ν + i A_ν,

where the A_ν constitute the coordinate expression of the local U(1) gauge connection. Ω_ν is given by 

    Ω_ν(x) = -i/4ω_abν(x)σ^ab

where σ^bc = i[γ^b,γ^c]/2 and ω_abν are the spin connection coefficients given by 

    ω^a_b_ν = e^a_μ∂_ν (e^μ_b) + e^a_μe^σ_bΓ^μ_σν.

Note that ω_abν = -ω_baν.
For a pedagogical introduction the reader is referred to <cit.>.

Specifically, for the metric given in (<ref>), the non-vanishing Christoffel symbols for the metric  are readily calculated to be
Γ^0_00= Φ_t, Γ^0_10= Φ_x Γ^0_11= Ψ_te^2(Ψ-Φ), Γ^0_22=Ψ_te^2(Ψ-Φ), Γ^1_00=Φ_x e^2(Φ-Ψ), Γ^1_10=Ψ_x, Γ^1_11=Ψ_x, Γ^1_12=Ψ_y, Γ^1_22=- Ψ_x, Γ^2_00=- Φ_y e^2(Φ - Ψ), Γ^2_11=-Ψ_y, Γ^2_21=Ψ_x, Γ^2_22=Ψ_y. Here we denoted f_α :=∂_α f and omitted symbols that may be recovered from the given ones by index symmetry.
The vielbein e^a_μ is readily calculated to be  e^0_0 = e^Φ e^1_1= e^2_2= e^Ψ. 
Since the metric is diagonal, we have for the non-vanishing spin connection components that numerically
ω^a_b_c =  Γ^a_bc if a ≠ b.
From these, one obtains

    Ω_0    = 1/4(ω_010 [γ^0, γ^1] + ω_020 [γ^0, γ^2] ) 
       =1/4(Φ_xe^Φ-Ψ [γ^0, γ^1] + Φ_ye^Φ-Ψ [γ^0, γ^2] ), 
    Ω_1    = 1/4(ω_011 [γ^0, γ^1] + ω_211 [γ^2, γ^1] ) 
       =1/4(Ψ_te^Ψ-Φ [γ^0, γ^1] + Ψ_y [γ^2, γ^1] ), 
    Ω_2    = 1/4(ω_022 [γ^0, γ^2] + ω_122 [γ^1, γ^2] )
       =1/4(Ψ_te^Ψ-Φ [γ^0, γ^1] + Ψ_x [γ^1, γ^2] ).


Working with local gamma matrices {γ^a}, we may write down the Dirac equation as 

    ie^-Φ∂_tψ -γ^0 mψ  + ie^-Ψγ^0γ^1∂_x ψ  + ie^-Ψγ^0γ^2∂_y ψ
       +i(e^-Ψγ^0γ^1[γ^0,γ^1]/4+ e^-Ψγ^0γ^2[γ^0,γ^2]/4)Ψ_te^Ψ-Φψ
       +i( e^-Φ[γ^0,γ^1]/4Φ_x e^Φ-Ψ + e^-Ψγ^0 γ^2[γ^1,γ^2]/4Ψ_x  )ψ
       +i( e^-Φ[γ^0,γ^2]/4Φ_y e^Φ-Ψ + e^-Ψγ^0 γ^1[γ^2,γ^1]/4Ψ_y  )ψ  =0

Choosing γ^0 = σ_z, γ^1=  - iσ_x and γ^2=  iσ_y, this equation becomes

    i(∂_t+Ψ_t)ψ= -iσ_ye^Φ-Ψ(∂_x + 1/2Ψ_x +1/2Φ_x )ψ
       -iσ_xe^Φ-Ψ(∂_y + 1/2Ψ_y +1/2Φ_y )ψ  +σ_ze^Φ m ψ.


Lastly, we remark that while it is well known that under a general diffeomorphism, a Dirac spinor can be taken to transform as a collection of  scalars <cit.>, if we work in coordinates and apply  a change of coordinates from {x^μ} to {y^μ}, then not only are we applying a diffeomorphism but we are also implicitly changing the frame in each tangent space with the help of which we express our Dirac spinor. This amounts to going from the basis {∂/∂ x^μ} to the basis {∂/∂ y^μ}. In case the change can be represented by a proper Lorentz transformation in each tangent space, standard arguments yield that this can be lifted to the spinor representation as

    ψ(x) ⟼Λ_1/2(ω_ab(Λ(x)))Ψ(Λ^-1(x)x),

where Λ is the corresponding matrix of the fundamental representation,

    Λ_1/2(ω_ab(Λ(x))) := exp( i/2ω_ab(x)i/4[γ^a,γ^b]).

Λ_1/2(·) denotes the spin 1/2 representation of the group element Λ and the antisymmetric ω_ab characterizes our Lorentz transformation  <cit.>.




 §.§ Effective coupling constants: Justification for neglecting fast oscillating contributions

Our starting point is (<ref>),

    k_[(a-1)↔ a]b^eff :=1/T∫_z_0^z_0 +T k_[(a-1)↔ a]b(z)dz,

where T:= 2π/ω, repeated here for convenience.
We assume that the scale of variation in the z-direction of the modes c_ab(z) is much larger than 1/α_j for all j. 
Furthermore we chose

    α_j=q_j ω,

with q_j positive natural numbers. 
We  expand the integrand using the generating function of the Bessel functions of the first kind

    exp(i σsin(x) ) = ∑_n∈Z J_n(σ)e^inx,

for each index j. 
To simplify notation, we write A^j:=A^j_(a-1)b-A^j_ab.
Interchanging summation and integration and using that 

    1/2 π∫_0^2 πexp(i∑_jq_jn_j t )dt = 
    1,    if ∑_jq_jn_j=0 
    
    0,    otherwise,

we see that

    k_[(n-1)↔ n]m^eff   = k_0  ∏_j=1^Λ J_0(  A^j/2 α_j)
       +k_0∑_(∑_j=1^Λ n_jq_j=0)
    not all n_j=0∏_j=1^Λ J_n_j(  A^j/2 α_j)^=:C_Λ.

Since the convergence in (<ref>) is absolute, we know that for fixed x

    ∑_n ∈Z |J_n(x)| < ∞.

Furthermore we note that |J_-n(x)|=|J_n(x)|. 
We now show by induction on Λ that by appropriately choosing the modulation speeds α_j, C_Λ can be made arbitrarily small.
For this, we define

    C_Λ^|·|:= ∑_(∑_j=1^Λ n_jq_j=0)
    not all n_j=0|∏_j=1^Λ J_n_j(  A^j/2 α_j)|.

Clearly
|C_Λ| ≤ C_Λ^|·|
and it suffices to show that  C_Λ^|·| can be made arbitrarily small.
For this choose the natural numbers in (<ref>) as

    q_j := 2^jΓ,

with a natural number Γ.
Since the case Λ=1 is trivial consider the case Λ=2 as the base case. There we have

    C_2^|·|=∑_n_2≠ 0|J_-2^Γ n_2( A^1/2 α_1) J_ n_2( A^2/2 α_2) |.

Thus it follows that

    C_2^|·|   ≤(∑_n≠ 0|  J_ n( A^2/2 α_2)  |   ) 
       ·(∑_k≠ 0|  J_ (k 2^Γ )( A^1/2 α_1)  |   ) Γ→∞⟶ 0.

Hence we may make C_2 arbitrarily small by choosing Γ large enough.

Now assume that C_Λ^|·|Γ→∞⟶ 0 for Λ≤ p ∈N and consider the case Λ = p+1.
To proceed, note that one may split the sum in C_Λ^|·| into a part where no n_j =0 and a part where always at least one n_j=0. The latter is then a sum of Λ =p+1 terms of the form  |J_0(...)|  C_p^|·| with some of the summands in C_p^|·| set to zero. These terms have vanishing magnitudes by assumption and the remaining term is 

    ∑_(∑_j=1^p+1 n_j2^jΓ=0)
     n_1,...,n_p+1≠ 0|∏_j=1^p+1 J_n_j(  A^j/2 α_j)| =: C̃^| ·|.

We note that summands where we have  |n_j| ≤ 2^Γ - 1 for all 1 ≤ j ≤ p+1 do not occur, since in this case we have

    |∑_j=1^pn_j2^jΓ|    ≤ (2^Γ-1)∑_j=1^p2^jΓ
       =(2^Γ-1)(2^Γ(p+1)-1/2^Γ - 1)
       = 2^Γ(p+1)-1 < 2^Γ(p+1)
       < |n_p+1|2^Γ(p+1),

violating the summation condition

    ∑_j=1^p+1 n_j2^jΓ=0.

Thus we in fact have

    C̃^| ·|≡∑_(∑_j=1^p+1 n_j2^jΓ=0)
    some |n_j| > (2^Γ - 1)
     n_1,...,n_p+1≠ 0
    |∏_j=1^p+1 J_n_j(  A^j/2 α_j)|.

From this we infer

    C̃^| ·|≤∑_k=1^p+1∑_(∑_j=1^p+1 n_j2^jΓ=0)
     |n_k| > (2^Γ - 1)
     n_1,...,n_p+1≠ 0
    |∏_j=1^p+1 J_n_j(  A^j/2 α_j)|.

But then clearly (implicitly keeping the summation over all n_j  in the first row but not in the the second and third)

    C̃^| ·|≤   ∑_k=1^p+1∑_ |n_k| > (2^Γ - 1)|∏_j=1^p+1 J_n_j(  A^j/2 α_j)|
    
    =   ∑_k=1^p+1( ∑_ |n_k| > (2^Γ - 1)|J_n_k(  A^k/2 α_k)|)
       ×∏_j≠ k( ∑_ n_j|J_n_j(  A^j/2 α_j)|).

This tends to zero as Γ tends to infinity.




 §.§ Vanishing diagonal couplings

In this appendix, we explain a possible method to realize vanishing diagonal couplings in a two dimensional square lattice. To avoid a cluttered notation, we will mainly resort to a pictorial representation. 
Our starting point are the approximate effective coupling constants (<ref>) and (<ref>) repeated here for convenience:

    k_[(a-1)↔ a]b   ≈ k_0  ∏_j=1^Λ J_0(  A^j_(a-1)b-A^j_ab/2 α^j)
    
    k_a[(b-1)↔ b]   ≈ k_0  ∏_j=1^Λ J_0(  A^j_a(b-1)-A^j_ab/2 α^j)


We set Λ =3 and define χ:= A/2 α^1 = B/2 α^2 such that J_0(χ)=0; E.g. χ≈ 2.405.
We then assign amplitudes according to

    [
    -A         -1/2A        0     -   1/2A   
                             
    1/2 A         A        1/2A        0    
    ],
 
for the first modulation and 

    [1/2 B         0        1/2B        0    
                             
    
    0         -1/2B        -B        -1/2B   
    ],
 
for the second modulation, periodically extended through the lattice.
As is easily verified directly this implies that for the coupling between sites that are diagonally adjacent, there is always at least one Bessel function yielding the value zero.
The coupling in y-direction is so far then given by 

    k_y = k_0 J_0(3/2χ) J_0(1/2χ).

The coupling in x direction is given by 

    k_x=k_0 (J_0(1/2χ))^2.

Note that k_y<0<k_x for ξ≈ 2.405.

To rectify this, we can use the third modulation.
We choose C/2α_3 such that 

    J_0(C/2α_3) < 0,

 and prescribe amplitudes according to

    [
    C         C        C        C    
                             
    
    0         0        0        0  
    ].
 

This then leaves k_x unchanged and ensures 0<k_x,k_y.
The remaining mismatch between the effective couplings can be absorbed into  different discretisation lengths in the two directions.


00

Aharonov Y. Aharonov, and D. Bohm, .  Phys. Rev. 115  485   (1959)
DeserS. Deser, R. Jackiw and  G. 't Hooft,  Ann. Phys. 152, 220 (1984)




Bezerra VB. Bezerra  Phys. Rev. D  35, 2031 (1987).


Ford L. H. Ford and   A. Vilenkin,  J. Phys. A 14, 2353 (1981).


Vilenkin0 A. Vilenkin,   Phys. Rev. D 23, 852 (1981). 


Burges  C. J. C. Burges, Phys. Rev. D 32, 504 (1985). 



Peschel U. Peschel, T. Pertsch, and F. Lederer, Opt. Lett. 23, 1701 (1998).
Morandotti R. Morandotti, U. Peschel, J. S. Aitchison, H. S. Eisenberg, and Y. Silberberg, Phys. Rev. Lett. 83, 4756 (1999).
Lahini Y. Lahini et al., Phys. Rev. Lett. 100, 013906 (2008).
Longhi11 S. Longhi, Appl. Phys. B 104 453 (2011).
Longhi12 S. Longhi and G. D. Valle, New J. Phys. 14, 053026 (2012).
Crespi A. Crespi, G. Corrielli, G. Valle, G. D. Valle, R. Osellame, and S. Longhi, New J. Phys. 15, 013012 (2013).
Keil R. Keil et al., Nat. Commun. 4, 1368 (2013).
Rodriguez13 B. M. Rodríguez-Lara, F. Soto-Eguibar, A. Zárate Cárdenas, and H. M. Moya-Cessa, Opt. Express 21, 12888 (2013).
Rodriguez14 B. M. Rodríguez-Lara and H. M. Moya-Cessa, Phys. Rev. A 89, 015803 (2014).
LeeAngelakis14 C. Lee, A. Rai, C. Noh, and D. G. Angelakis, Phys. Rev. A 89, 023823 (2014).
Marini14 A. Marini, S. Longhi, and F. Biancalana, Phys. Rev. Lett. 113, 150401 2014).
RaiAngelakis15 A. Rai, C. Lee, C. Noh, and D. G. Angelakis, Sci. Rep. 5, 8438 (2015).
Keil15 R. Keil et al., Optica 2, 454 (2015).

Longhi10a S. Longhi, Opt. Lett. 35, 235 (2010).
Longhi10b S. Longhi, Phys. Rev. B 81, 075102 (2010).



Dreisow10 F. Dreisow, M. Heinrich, R. Keil, A. Tünnermann, S. Nolte, S. Longhi and A. Szamzeit, Phys. Rev. Lett. 105, 143902 (2010).



Dreisow12 F. Dreisow, R. Keil, A. Tünnermann, S. Nolte, S. Longhi and A. Szamzeit, Euro. Phys. Lett. 97, 10008 (2012).

                       

Setare M. R. Setare, P. Majari, C. Noh and Sh. Dehdashti, J. Mod. Opts. 66, 1663 (2019).


Boada O. Boada, A Celi, J I Latorre and M Lewenstein, New J. Phys 13, 035002 (2011).



Kosior A. Kosior, M. Lewenstein and A. Celi,  SciPost Phys. 5, 061 (2018).

Celi  A. Celi, Eur. Phys. J. ST 226, 2729 (2017).




Koke
C. Koke, C. Noh and D. G. Angelakis,  Ann. Phys. 375, 162 (2016).




Efremidis N. K. Efremidis, P. Zhang, Z. Chen, D. N. Christodoulides, C. E. Ruter, and D. Kip, Phys. Rev. A 81, 053817 (2010)

Szameit09 A. Szameit et al. Phys. Rev. Lett. 102, 153901 (2009).  

Birell N.D. Birell and P.C.W. Davies, Quantum Fields in curved Spacetime, Cambridge University Press  (Cambridge, 1982)


Szameit A. Szameit and S. Nolte, J. Phy. B: At. Mol. Opt. Phys. 43, 163001 (2010).






VilenkinA. Vilenkin, E. P. S. Shellard, Cosmic strings and other topological defects; Cambridge University Press (Cambridge, 2000)



Thaller B. Thaller Visualizing the kinematics of Dirac wave packets,  quant-ph/0409079 (unpublished)  (2004).

Lienert M. Lienert, Wave equations of relativistic quatnum mechanics,   Eberhard-Karls Universität Tübingen, (unpublished) (2018).

	
Peskin M.E. Peskin,  and D. V. Schroeder, An Introduction to Quantum Field Theory;  Addison-Wesley (Reading (USA), 1995)


Nakahara M. Nakahara, Geometry, topology and physics;
 Taylor & Francis (Boca Raton (USA), 2003).



Carlip S. Carlip,   J. Korean Phys. Soc. 28, 447 (1995)  

Garcia
 A. Garcia, F. W. Hehl, C. Heinicke , A. Macias,  Class. Quant. Grav. 21, 1099  (2004).


Sitenko2010  Y. Sitenko  J. Phys. A43, 354014  (2010).

Isabelstring I. Fernández-Núñez, O. Bulashenko Phys. Lett. A 380, 2897 (2016) 


Zwiebach B. Zwiebach, A First Course in String Theory,  Cambridge University Press  (Cambridge, 2009)









