pISSN 0374-4914 eISSN 2289-0041

## Research Paper

New Phys.: Sae Mulli 2021; 71: 1096-1104

Published online December 31, 2021 https://doi.org/10.3938/NPSM.71.1096

## One-dimensional Variable-mass Dirac Equation and Spinor slow Llght

Changsuk NOH*

Department of Physics, Kyungpook National University, Daegu 41566, Korea

Correspondence to:cnoh@knu.ac.kr

Received: September 24, 2021; Accepted: October 12, 2021

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License(http://creativecommons.org/licenses/by-nc/3.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

The one-dimensional variable-mass Dirac equation is connected to various models used throughout many branches of physics. An analog simulation of the equation in a spinor slow light system allows experimental realizations of such models. This work concentrates on an interesting model of historical importance that led to a prediction of charge fractionalization, which in turn occurs due to the presence of a topologically protected zero-energy mode. After describing how the model can be realized in a spinor slow light system, the current work explains how the presence of the zero-energy mode can be verified from the dynamics of the spinor slow light.

Keywords: Dirac equation, Quantum optics, Dark state polaritons

The Dirac equation was the first successful attempt at merging quantum mechanics and relativity [1]. The relativistic wave equation predicted paradoxical effects such as Zitterbewegung, a rapid oscillatory motion of a Dirac particle, and Klein’s paradox, tunneling of a Dirac particle through an infinitely long barrier. The cause of these effects lies in the fact that the Dirac equation allows both positive and negative energy states, or particle/ antiparticle pairs as discovered by Dirac himself. The subsequent development of quantum field theory explained these paradoxical effects, but the equation itself and its predictions continued to fascinate many generations of physicists.

An analog simulation, or an emulation, of the Dirac equation was first proposed in a trapped ion system in 2007 [2]. The one-dimensional Dirac equation was implemented experimentally soon after, and both Zitterbewegung [3] and Klein’s paradox [4,5] have been observed. This has sparked an interest in emulations of the Dirac equation in other systems such as coupled waveguides [69] and dark state polaritons [10,11]. These were followed by many proposals and experiments on emulations of relativistic models and effects, which are too numerous to mention here (see [1217] for a selection).

Usefulness of the Dirac equation is not limited to the description of a fermionic particle. By allowing the mass term to vary in space (or even in time) the model’s applicability extends to many other areas of physics including condensed matter physics and nuclear physics. Pointing out such a connection is one of the aims of this work. For simplicity and for the ease of experimental implementation, we stick to the one-dimensional Dirac equation, which, despite the limitation, gives rise to plenty of interesting physics as we will see. Only one-dimensional models will be considered in the rest this work and the Dirac equation will always be in one-dimension.

A convenient starting point for a derivation of the onedimensional Dirac equation is the covariant form

$iγμ∂μψ=mcψ,$

where {γµ, γν} = 2ηµν with ηµν the metric in flat spacetime. The latter is the two-by-two Pauli matrix σz in one spatial dimension and as a consequence, other Pauli matrices can be used to represent the gamma matrices γµ. This implies that the spinor has two components instead of the four needed in three-dimensional space. Using γ0 = σy and γ1 = x, one obtains

$i∂tψx,t=icσz∂xψx,t+mc2∂yψx,t,$

with a two-component spinor ψ(x, t) = (ψ1(x, t), ψ2(x, t))T. This equation will be called the variable-mass Dirac equation when the mass term mc2 is allowed to vary in space (and time).

This work provides brief descriptions of various models associated with the variable-mass Dirac equation and show how some of their predictions can be emulated using spinor slow light (SSL). Section II gives a gentle introduction to the so-called dark state polaritons (see Ref. [18] for a short pedagogical introduction), concentrating on how they can be coerced to form an SSL, which follows the dynamics of the Dirac equation. It also explains the preparation, evolution, and measurement stages of the emulation process and sets the stage for the rest of the work. Section III introduces the Jackiw- Rebbi model, whose implementation in an SSL system has already been proposed in [19]. Whereas the latter work concentrated on a transmission-and-reflection scenario, this work concentrates on a loading-and-releasing scenario. Section IV discusses three other models related to the variable mass Dirac equation, namely polyacetylene, the random mass Dirac model, and the Lorentz scalar potential.

### II. Spinor slow light

That the so-called dark state polaritons (DSPs) can be manipulated to simulate the physics of the Dirac equation has been noticed in Refs. [10,20] and experimentally realized in [21]. The basic setup consists of atoms coupled to an effectively one-dimensional waveguide as illustrated schematically in Fig. 1(a). Atomic levels have a double-tripod structure as shown in Fig. 1(b). In order to understand this setup, let us first consider a simpler atomic structure of the lambda-type configuration depicted in Fig. 2(a). When a probe field resonant with the atomic transition |g⟩ ↔ |e⟩ enters the waveguide, it is strongly absorbed by the medium.

Figure 1. (Color online) (a) Schematic representation of the spinor slow light setup and (b) associated double tripod atomic level structure.

Figure 2. (Color online) Atomic level schemes of the lambda type that realizes (a) the slow light and (b) the stationary light.

The situation changes completely when a strong control field, driving the transition |e⟩ ↔ |s⟩, is introduced. Due to a destructive interference, a phenomenon called electromagnetically induced transparency (EIT) sets in: The probe field passes freely through the medium and does so at an ultra-slow speed [22]. In the medium, a collective excitation of light and matter called a dark state polariton is formed [22], which follows the dynamics governed by the equation

$∂t+υg∂zψ^x,t=αψ^z,t.$

$ψ^$ is the dark state polariton operator, which is proportional to the probe field operator $E^$ under the condition of EIT, i.e., a small group velocity. The proportionality factor α depends on system parameters such as the atom-field coupling strength (g), atomic density (n), control field strength (Ω), and the two photon detuning (δ). The precise form of α is not needed for the purpose of this work. The group velocity is given by vgvΩ2/(Ω2 + 2πg2n), where v is the speed of light in an empty waveguide.

By introducing a counter-propagating control fields Ω1,2 to the above system as shown in Fig. 2(b), two counter propagating DSPs, $ψ^$1,2, are formed, which combine to form a single polariton field. The latter has zero group velocity and is called the stationary light [23], which can be formed as follows. First, a traveling pulse of light E1 is introduced to the system with the control field Ω1 turned on. Due to the EIT condition, the pulse travels slowly through the medium. When the pulse has completely entered the medium, the control field is turned off adiabatically, which reduces the group velocity to zero. The pulse is now stored as atomic coherences between the states |g⟩ and |s⟩ and can be retrieved at will by turning the control field back on. A DSP system can therefore be used as a quantum memory device as experimentally demonstrated in Ref. [24]. A longlived stationary light is produced when the control fields Ω1 and Ω2 are slowly turned on after the pulse has been loaded.

To obtain a spinor slow light, consisting of twocomponents, additional ground (or metastable) states are needed [10,20]. The double tripod setup proposed in [20] is shown in Fig. 1(b). Such a configuration yields two DSPs, Ψ1,2, that obey an effective equation

$i∂t+iυ0σz∂z−δσyΨz,t=0,$

where Ψ = (Ψ12)T , v0 is the group velocity, and δ is the two-photon detuning as depicted in the figure. Note that this takes the same form as Eq. (2) so the variablemass Dirac equation can be emulated by allowing the two-photon detuning δ to vary in the spatial dimension z or even in time t. Only spatially varying mass will be considered in this work, but time-varying models can be related to interesting physics in curved 1+1 dimensional spacetime as shown in Ref. [8].

### III. Solitons and the Dirac equation

In its essence, the Jackiw-Rebbi (JR) model is described by the Dirac equation with a specific type of position-dependent mass profile [25]. It was used by the authors to predict charge fractionalization, before fractional quantum Hall effect was discovered. Essentially the same effect has been found independently in studying the electron-phonon coupling in polyacetylene [26]. Before we delve into the details of the model, let us start with a brief remark on charge fractionalization, the details of which can be found in a review article [27].

The JR model arises as an approximation to the quantum field theory for an electron coupled to bosonic fields. The equations of motion for the bosonic fields allow soliton solutions when treated classically, which, in the approximate Dirac equation, acts as a position-dependent mass term with the soliton profile. A fractionalization of charge occurs in this model because it allows a zeroenergy mode, or a zero-mode for short. The ground state of the quantum field can either contain the zero-mode or not, and the charge conjugation symmetry dictates that these two degenerate ground states must have opposite charges. Now, because only one electron may occupy the zero-mode, the latter must have a fractional charge of e/2 or -e/2. Indeed, a detailed calculation of the formal charge operator confirms this prediction. Surprisingly, the existence of a zero-mode is guaranteed irrespective of the details of the soliton and depends only on a global property. That is, the existence of the zero-mode is topologically protected.

### 1. Jackiw-Rebbi model

As noted above, charge fractionalization is ultimately connected to the presence of the zero-mode, which is topologically protected by global properties of the soliton field. While the polariton emulator cannot exhibit fractionalization explicitly because the DSP is not a fermionic field, the existence of the zero-mode and its topological robustness can be probed as shown in Ref. [19]. This section explains the basic theory of the zero-mode and describes how it can be observed in the spinor slow light system. While the original proposal has concentrated on optical transmission measurements [19], this work concentrates on the dynamics of the spinor.

In the JR model, a Dirac particle is coupled to a real scalar field φ(z) as described by

$i∂tΨ=σzcpz+mc2κϕzσyΨ,$

where φ(z) obeys the Klein-Gordon equation

$∂μ∂μϕz+λ22κ2 κ2 −ϕz2 2=0.$

The ground state solutions are φ(z) = ±κ and the degeneracy implies a soliton solution that interpolates between −κ at z = −∞ and κ at z = ∞ (the corresponding anti-soliton also exists). That is, there exists a (soliton) solution of the form

which is localized around z = 0. Plugging this into Eq. (5) yields the zero-energy solution

where σzσyχ = −x = −iχ holds, yielding χ = [1, 1]T / $2$. That this is the zero-energy solution can be easily checked by substituting the solution into Eq. (5). The soliton solution and the corresponding zero-mode are depicted in Fig. 3(a).

Figure 3. (Color online) The spatial part of the zero-mode spinor, Ψ0(z) shown as the solid (blue) curves along with the ‘mass’ profile φs(z) shown as (black) dashed curves. (a) For the soliton solution φ(z) = κ tanh(λz) and (b) for the linear profile φ(z) = λz. Localization scale is on the order of λ = 1. Other parameters are mc = 1 and κ = 0.5.

Before we move on to the topic of emulation using SSL, a discussions on a few important aspects of the zero-energy solution is in order. Firstly, the zero-mode solution only depends on the value of mc, upon fixing the parameters λ = 1, which determines the width of the solution. To understand the second aspect, consider a generic mass profile φ(z) instead of the soliton profile φs(z). Let us assume a constant φ(z). Then the integral on the first line of Eq. (7) diverges as z → −∞, so the solution is not normalizable. Therefore, a zero-energy solution exists only if φ(z) goes from a negative value at z = −∞ to a positive value at z = ∞. Lastly, a zeromode exists as long as φ(z) has these limiting values. While the detailed shapes of the solution depends on the mass profile, the existence of the zero-mode depends only on the fact that φ(z) has the correct boundary condition. This is evident for a linear profile φ(z) = z, which is illustrated in Fig. 3(b) along with its zero-energy solution. Although such φ(z) is not a valid solution of the Klein- Gordon equation, it serves as a nice example to show the topological nature of the zero-energy solution.

### 2. Analog simulation with spinor slow light

Due to the similarity between Eqs. (4) and (5), the SSL system follows the dynamics of the Dirac equation with an effective speed of light v0 = c and an effective mass term δ = $mc2κ$φ(z). The JR model can be emulated if the two-photon detuning δ can be made to follow the soliton profile, which can be achieved by introducing an auxiliary field to shift the energy levels appropriately [21]. Having prepared such a system, how can one ‘observe’ the zero-mode? Two possible schemes have been proposed in [19].

In the first scheme, one utilizes the intrinsic nonequilibrium property of the quantum optical systems by driving the waveguide with a laser field [28] and observing the transmission spectrum [19]. Transmission spectrum shows a peak at zero energy, revealing the presence of the zero-mode.

The second scheme, which has been only briefly touched upon in [19], focuses on the time evolution of a spinor. To understand this scheme, let us assume that the zero-energy state has been prepared at time t = 0. Then because it is an eigenstate of the Dirac Hamiltonian, the state remains unchanged for t > 0. Because preparing the exact eigenstate will be difficult, let us consider an initial Gaussian state. That is, ΨGauss(z, 0) = N exp[−z2/2σ2]χ where N is a normalization factor. This class of states can be used to closely mimic the ideal zero-energy state as illustrated in Fig. 4. For λ = 1, σ = 1.2 yields the state overlap of ∫ dzψGauss(z, 0)ψ0(z, 0) ≈ 0.997, which is equivalent to the overlap between two quantum states ⟨ψ>12⟩. ψGauss and ψ0 are the spatial parts of the Gaussian state and the zero-energy state, respectively.

Figure 4. (Color online) Comparison between a Gaussian function with σ = 1.2 (solid blue curve) and the spatial part of the zero-mode state with λ = 1 (dashed red curve).

Under time evolution, the initial Gaussian state exhibits completely different behavior depending on whether the mass term has a soliton profile or not. In Fig. 5(a) the mass term δ has the soliton profile tanh(z), in which case the total intensities of each component of the spinor remain unchanged (changes are too small to be noticed). By total intensity, we mean the integral of the square of the absolute value of a spinorfcomponent, i.e., Ii ≡ ∫ dzi(z, t)|2. For a constant mass profile, say δ = 1, the intensities fluctuate as shown in Fig. 5(b). Therefore, the existence of the zero-mode can be inferred from the time evolution of the intensities, which can be observed in the SSL system by retrieving the light pulses.

Figure 5. (Color online) Total intensities of the spinor components as functions of time for an initial Gaussian spinor with σ = 1.2. Solid blue curves represent I1 while dotted orange curves represent I2. (a) When the mass term is given by the soliton mass profile δ(z) = tanh(z). (b) For a constant mass δ = 1.

Actually, this is not the whole story and Fig. 5(a) can be deceiving. For example, for σ = 0.5 the fidelity drops to ~ 0.85, but the figure remains the same. In fact, the total intensities do not fluctuate at all for an initial Gaussian state regardless of the value of σ. To observe differences in the dynamics, one needs to look at the intensity profiles as illustrated in Fig. 6. Panel (a) displays the density plots of the intensity profiles |Ψi(z, t)|2 for the initial Gaussian state with σ = 1.2. There are slight modulations in the profile but they are well confined to the central region. Upon reducing σ to 0.5, the fidelity drops and the two components scatter in the opposite directions as shown in panel (b). However, even in this case, a significant fraction of the spinor eventually gets trapped in the zero-mode as evidenced by the central peaks after time t ≳ 3,

Figure 6. (Color online) Time evolution of the spinor. Density plots of the intensities of the spinor components |Ψ1|2 and |Ψ2|2. The brighter the color, the higher the intensity. Initial Gaussian state with (a) σ = 1.2 and (b) σ = 0.5.

An experiment to verify the presence of the zero-mode can be run as follows. First, a Gaussian state is prepared initially. This can be achieved either by using the EIT technique described earlier to store two counter propagating light pulses as atomic coherences or by using Raman transitions to create the coherences directly. Second, all four control fields Ωij with i, j = 1, 2 are turned on, which effectively creates the required Dirac Hamiltonian. The JR model can be emulated if the frequencies of the control fields vary in such a way to make δ follow the soliton profile. Another method is to use auxiliary electric or magnetic fields to detune the states |si⟩ [21].

Third, after waiting for a certain period of time t, turn off the control fields and retrieve the light using the conventional EIT technique. Intensities of the retrieved fields correspond to the intensities of the spinor components. Therefore, by measuring how the intensities of the retrieved fields change in time, the intensity distributions |Ψi(z)|2 for a given run time t can be measured. By repeating the experiment for different values of t, an equivalent of Fig. 6 can be constructed.

There can be many complications in an actual experiment, one of which is that the total intensity decays in time [21,29]. This occurs for many reasons, such as couplings to other states in an actual atom and instability of the states |s1⟩ and |s2⟩. Generally, the two spinor components have different decay rates and by assuming that the decay rates are uniform in z, the equation of motion becomes

$i∂t+iυ0σz∂z−δσy−iΓΨz,t=0,$

where Γ is a 2 by 2 matrix with diagonal components γ1 and γ2, corresponding to the decay rates of Ψ1 and Ψ2, respectively. This causes the total intensities of the individual components to decay, limiting the evolution time, but the effects of the zero-mode can still be observed from the initial time evolution. This is illustrated in Fig. 7, where an asymmetrical choice γ1 = 0.1 and γ = 0.2 was made. Only |Ψ1(z, t)|2s are shown, for an initial state that is very close to (left) and not quite the same as (right) the zero-mode. |Ψ2(z, t)|2 decays a little faster. A comparison to Fig. 6 clearly demonstrates the effect of the dissipation.

Figure 7. (Color online) Time evolution of |Ψ1(z, t)|2 in the presence of loss. The loss parameters are γ1 = 0.1 and γ2 = 0.2. (a) σ = 1.2 and (b) σ = 0.5.

Lastly, the topological character of the zero-mode can be verified by varying the details of the spatial mass profile. For example, upon changing the value of λ, one can always find a zero-energy solution via Eq. 7, which can be verified by a judicial choice of σ. In an experiment, it is actually easier to choose the linear mass profile considered in Sec. III 1, because the zero-mode is a Gaussian in such a case [21].

### IV. Dirac equation in other areas of physics

The variable-mass Dirac model is not only useful in high-energy physics, but also in other branches of physics such as condensed matter and nuclear physics. The purpose of this section is to point out three such cases: Polyacetylene, random mass Dirac model, and Lorenz scalar potential. We limit ourselves to the most basic aspects of each topic and briefly discuss how the relevant physics can be emulated in a SSL setup.

### 1. Polyacetylene

Polycetylene can be modeled as electrons confined to move in a one-dimensional chain [27]. Phonon modes of the atoms are coupled to the electrons because the hopping strength depends on the distance between the sites. In such a system, the configuration of equally spaced sites does not support the ground state. Instead, it supports two degenerate ground state configurations with broken reflection symmetry. This is called the Peierls instability [30]. A soliton is formed when the two configurations coexist in a single chain, at the location at which the configuration changes from one to the other. A localized electron state forms around the site, whose energy lies between the gap created by the electron-phonon coupling.

The situation is analogous to the Jackiw-Rebbi model, where a zero-mode forms between the mass gap of the Dirac electrons. This can be established formally by taking the continuum limit of the electron hopping model and using the Jordan-Wigner transformation, which turns out to yield a slightly more generalized version of Eq. (5). Therefore, analogous physics to the JR model occur in polyacetylene, which can be observed in the same SSL setup described above. One complication that arises in polyacetylene is the spin of the electron, which is absent in the one-dimensional Dirac particle. In short, the soliton is associated with a charge e and spin 0, so the singly-occupied state has charge 0 and spin ±1/2 [27]. This is an anomalous spin-charge relation arising due to charge fractionalization.

### 2. Random mass Dirac model

Random mass Dirac model is a variant of the Jackiw- Rebbi model, where solitons are assumed to take the form of a step-function (a kink) and are randomly distributed in space. The mass term in the Dirac equation thus randomly flips its sign. This is called the telegraph signal and is described by

$ϕx=∏isgnx−x i,$

where sgn denotes the sign function and xi denotes the randomly distributed locations of the kinks. Interest in this model arose from the realization that the Dirac equation describes the low-energy physics in the spin sector of the two chain spin-ladder [31] and a spin-Peierls system [32].

The spin-Peierls model describes an antiferromagnetic spin-1/2 chain with an additional dimerization term [32]:

$HSP=∑nJ+−1 nΔ nS→ nS→n+1.$

Δn determines the strength of the phonon-mediated dimerization and is related to the mass term in the continuum limit. When impurities are introduced, domain walls are formed at their locations, affecting Δn, which in turn introduces a kink in the mass term. The lowenergy physics of this model is described by the random mass Dirac model.

One of the characteristics of the model is a long-range correlation in the mid-gap (zero-mode) state [32]. Because the localization length of a state is proportional to ln(1/ϵ), where ϵ is the energy of the state, the mid-gap state has a diverging localization length. This in turn leads to a polynomial decay in the impurity-averaged intensity-intensity correlation function

$IxI0¯∝x−3/2,$

at length scales greater than the mean free path.

An emulation in terms of our polaritonic setup starts from noting that the Gaussian distribution of the mass term, instead of the telegraph signal described earlier, exhibits the same scaling in the correlation function [32]. Note that the mass term does not follow a Gaussian profile in space, but has a random value picked out from a Gaussian distribution. In an emulation, the two-photon detuning should be varied randomly in space according to a Gaussian distribution. After loading and evolving an initial state following the same procedure introduced in the previous section, the intensity correlation function I(x)I(0) can be obtained by measuring the intensity correlations in the spinor components. The latter converts to time-correlations in the retrieved field, that is, |Ψi(t)|2i(0)|2 where i = 1, 2. Finally, averaging over many impurity realizations yields the desired correlation function.

We note that the random mass Dirac model can be realized in the SSL system was already noted by Unanyan et al. [10] and a related model with a random vector potential has been studied in the setting of cold atoms [33]. Even an experimental demonstration was performed in a waveguide arrays setup [7].

### 3. Lorentz scalar potential

As another application of the variable-mass Dirac equation, imagine a Dirac field coupled to a Lorentzscalar potential. Such a case was first studied as a spinoff of the MIT bag model, devised as a phenomenological model for quark confinement [34, 35]. Yet another area that uses similar models is nuclear physics (see e.g. [36]). The Hamiltonian describing a Dirac field interacting with both a scalar potential Vs(x) and a vector potential Vv(x) can be written as

$i∂tΨ=αpz+Vυx+m+VsxβΨ.$

The vector potential is known to exhibit Klein’s paradox when it is strong enough to create a particle and an antiparticle [1]. This type of behaviour occurs because the positive and negative energy states see different potentials and does not occur for a scalar potential.

In the phenomenological study of quarkonium, an equal mixture of a vector and a scalar confining potential is used to obtain a best fit for the spin-orbit splitting [37]. One widely studied potential is a linearly increasing potential, V (x) = g|x|. A similar model is also used to study pseudospin doublets in nuclei [38].

The vector potential term can be introduced by slightly modifying the double tripod scheme in Fig. 1(b). Instead of choosing the equal and opposite two photon detunings as shown in the figure, choose different detunings δ1 and δ2. Then, following the analysis in Ref. [20], one obtains an extra term δ0Ψ(z, t) on the left hand side of Eq. (4), where δ1 = δ0+δ and δ2 = δ0δ. This allows an experimental realization of Eq. (12) with various combinations of vector and scalar potentials. In particular, fixing the vector potential term and tuning the strength of the scalar potential, by changing the two-photon detuning, would yield a crossover from the regime of the Klein paradox to confinement. Methods described in the previous section can be used to check such a transition experimentally, which will be clearly imprinted on the spatial profile of the SSL and thus the temporal profile of the retrieved field. Another option is to employ a continuous driving scenario, in which the scalar potential term can be changed continuously in time while the system is continuously driven, which will affect, in real time, the statistics of the output field.

In conclusion, a connection between the onedimensional Dirac equation and a spinor slow light system has been pointed out. The latter system was shown to be able to emulate a variable-mass Dirac equation which in turn is related to models in diverse branches of physics such as high-energy, condensed matter, and nuclear physics. In particular, the SSL setup realizing the Jackiw-Rebbi model was considered in detail, along with a possible experimental procedure to observed the presence of the topologically protected zero-mode.

This research was supported by Kyungpook National University Research Fund, 2019.

1. B. Thaller, The Dirac Equation. (Springer-Verlag, Berlin Heidelberg, 1992).
2. L. Lamata, J. León, T. Schätz and E. Solano, Phys. Rev. Lett. 98, 253005 (2007).
3. R. Gerritsma, G. Kirchmair, F. Zähringer, E. Solano, R. Blatt and C. F. Roos, Nature 463, 68 (2010).
4. R. Gerritsma, B.P. Lanyon, G. Kirchmair, F. Zähringer and C. Hempel et al, Phys. Rev. Lett. 106, 060503 (2011).
5. J. Casanova, J.J. García-Ripoll, R. Gerritsma, C.F. Roos and E. Solano, Phys. Rev. A 82, 020101(R) (2010).
6. F. Dreisow, M. Heinrich, R. Keil, A. Tünnermann and S. Nolte et al, Phys. Rev. Lett. 105, 143902 (2010).
7. R. Keil, J.M. Zeuner, F. Dreisow, M. Heinrich and A. Tünnermann et al, Nature Comm. 4, 1368 (2012).
8. C. Koke, C. Noh and D. G. Angelakis, Ann. Phys. 374, 162 (2016).
9. C. Koke, C. Noh and D. G. Angelakis, Phys. Rev. A 102, 013514 (2020).
10. R.G. Unanyan, J. Otterbach, M. Fleischhauer, J. Ruseckas, V. Kudriass̆ov and G. Juzeliūnas, Phys. Rev. Lett. 105, 173603 (2010).
11. J. Otterbach, R. G. Unanyan and M. Fleischhauer, Phys. Rev. Lett. 102, 063602 (2009).
12. O. Boada, A. Celi, J. I. Latorre and M. Lewenstein, New J. Phys. 13, 035002 (2011).
13. C. Noh, B. M. Rodríguez-Lara and D. G. Angelakis, New J. Phys. 14, 033028 (2012).
14. B.M. Rodriguez-Lara and H. M. Moya-Cessa, Phys. Rev. A 89, 015803 (2014).
15. R. Keil et al, Optica 2, 454 (2015).
16. A. Celi, Eur. Phys. J. Spec. Top. 226, 2729 (2017).
17. M.R. Setare, P. Majari, C. Noh and Sh. Dehdashti, J. Mod. Opt. 66, 1663 (2019).
18. C. Noh and D. G. Angelakis, Rep. Prog. Phys. 80, 016401 (2017).
19. D.G. Angelakis, P. Das and C. Noh, Sci. Rep. 4, 6110 (2014).
20. J. Ruseckas et al, Phys. Rev. A 83, 063811 (2011).
21. M. Namazi et al, arXiv:1711.09346 (2017).
22. M. Fleischhauer, A. Imamoglu and J.P. Marangos, Rev. Mod. Phys. 77, 633 (2005).
23. M. Bajcsy, A.S. Zibrov and M.D. Lukin, Nature 426, 638-641 (2003).
24. D.F. Phillips et al, Phys. Rev. Lett. 86, 783 (2001).
25. R. Jackiw and C. Rebbi, Phys, Rev. D 13, 3398 (1976).
26. W.P. Su, J.R. Shrieffer and A.J. Heeger, Phys, Rev. B 22, 2099 (1980).
27. R. Jackiw and J.R. Schrieffer, Nucl. Phys. B 190, 253 (1981).
28. P. Das, C. Noh and D.G. Angelakis, EPL 103, 34001 (2013).
29. M.-J. Lee, J. Ruseckas, C.-Y. Lee, V. Kudrias̆ov and K.-F. Chang, Nat. Commun. 5, 5542 (2014).
30. R.F. Peierls, Quantum theory of solids. (Clarendon Press, Oxford, 1955).
31. D.G. Shelton, A.A. Nersesyan and A.M. Tsvelik, Phys. Rev. B 53, 8521 (1996).
32. M. Steiner, M. Fabrizio and A.O. Gogolin, Phys. Rev. B 57, 8290 (1998).
33. S. Zhu, D. Zhang and Z.D. Wang, Phys. Rev. Lett. 102, 210403 (2009).
34. C.L. Critchfield, Phys. Rev. D 12, 923 (1975).
35. A. Chodos et al, Phys. Rev. D 9, 3741 (1974).
36. B.H.J. McKellar and G.J. Stephenson Jr., Phys. Rev. C 35, 2262 (1987).
37. P.M. Fishbane, S.G. Gasiorowicz, D.C. Johannsen and P. Kaus, Phys. Rev. D 27, 2433 (1983).
38. J.N. Ginocchio, Phys. Rep. 315, 231 (1999).