Subradiant emission from regular atomic arrays: universal scaling of decay rates from the generalized Bloch theorem
UUniversal Scaling of Decay Rates of Emitter Arrays in the Subradiant States
Yu-Xiang Zhang ∗ Niels Bohr Institutet, Københavns Universitet, 2100 København Ø, Denmark (Dated: June 12, 2020)Sub-wavelength arrays of emitters support subradiant excitations with decay rates scaled as N − α ,with N the number of emitters. In this Letter, we study one-dimensional arrays and find α = s + 1for the subradiant states with wave numbers close to an extremum ( k ex ) of the dispersion relation ω k , where the s -index is the order of the leading term of the expansion ω k − ω k ex ∝ ( k − k ex ) s . Therelation is universal regardless of the specific Hamiltonians. It also implies a new effect in systemswith multiple bands: the occurrence of band gap closing (e.g., at topological transitions), changesthe s -index hence surges the radiation from the subradiant states. This effect is demonstrated indimerized emitter arrays coupled to free space modes or an ideal waveguide. Subradiance is the phenomenon that radiation froman emitter ensemble is collectively prohibited [1]. Theopposite effect, superradiance, is well known for the N times enhancement of radiation rate ever since the Dickemodel [2], with N the number of emitters. However,the parallel problem that how the subradiant suppres-sion scales in N is not yet well understood. Recently, inone-dimensional (1D) arrays of equally spaced emitterscoupled to various light fields, subradiant states with de-cay rates scaled as N − are found [3–8] and consideredas a universal feature [4]. However, scalings of N − α with α > α , and to what extent they are truly universal.Answers to the questions will improve the understandingof the collective dynamics of emitter arrays [6, 10–21],and facilitate the engineering of subradiant states for po-tential applications in quantum memories [22, 23], energytransfer [24–26] topological photonics [27, 28], etc.To address them, consider a wave packet propagatingin a 1D emitter array coupled to some electromagneticfield. The duration for it to be trapped in the array isstill limited by the group velocity, even if it does notdirectly radiate. The group velocity will be slow if thewave number of the excitation is close to an extremum k ex of the dispersion relation ω k . Suppose that the orderof the extremum is s , i.e., we have the expansion ω k − ω k ex ∝ ( k − k ex ) s for k ≈ k ex . A larger s implies a slowergroup velocity (also, a flatter band [29]). Thus a relationbetween α and s is expected, while it is nontrivial tosubstantiate the universality.In this Letter, we shall show α = s + 1 for the single-excitation states of 1D arrays. This relation is estab-lished through a solvable toy model, which is then linkedto real systems by the primary perturbation formalism.The universality will be manifested from the approach.This relation implies a new effect: a change of the s -index(caused by, e.g., band gap closing at topological transi-tions) will induce a surge of the radiation of the sub-radiant states. We demonstrate this effect in dimerizedemitter chains in the 3D free α from 3 to 1 (numerically)or coupled to an ideal 1D waveguide (analytically), where a jump of α from 3 to 1 occurs at the critical point of theSu-Schrieffer-Heeger (SSH) type topological transition. Models.
We shall take the emitters arrays in the 3Dfree space as the exemplary system of this work. Thereof,the Born-Markov approximation works well hence thelight field can be traced out. It results in a non-HermitianHamiltonian that describes the effective emitter-emittercoupling induced by light [30]: H eff = − µ ω N X i,j =1 d ∗ i G ( x i − x j , ω ) d j σ † i σ j , (1)where µ is the vacuum permeability, ω is the transi-tion frequency between the emitter ground state | g i andexcited state | e i , d i and x i are the transition dipole mo-ment and coordinate of the i th emitter, σ i = | g i i h e i | , and G is the dyadic Green’s tensor of the Maxwell’s equa-tions. Dipole-dipole interactions induced by other lightfields are obtained by replacing G with the specifiedones [30]. As depicted in Fig. 1(a), we suppose the emit-ters are equally separated by d and polarized transverseto the array. We denote the system by DD-T in short-hand. DD-T was noticed to be more irregular because ofits 1 /r long-range interaction [4].For infinite chains, the Bloch’s theorem reveals a greatuniversality that the eigenstates are always the Blochstate, which is, if written for finite N , | k i = 1 √ N N X m =1 e ikx m σ † m | g g · · · g N i (2)where k ∈ [ − π/d, π/d ] (the first Brillouin Zone). Thespecified Hamiltonian determines only the eigenvalues.For the non-Hermitian H eff , the eigenvalue is written as ω k − iγ k /
2, where γ k is the decay rate of | k i (with infinite N ). In Fig. 1(b), ω k and γ k of DD-T with the resonantwave number k = 0 . π/d are depicted. It is notable that γ k = 0 if | k | > k , because no resonant light modes invacuum have so large wave numbers [4].Eigenstates of finite chains would have amplitude pro-files that look like standing waves, rather than the Blochstates. Nevertheless, their phase profiles help to associate a r X i v : . [ qu a n t - ph ] J un N −3 N −5 ( a ) x d ( b ) ω k γ k ( d )( c ) kd Figure 1. (a) 1D emitter array in 3D free space, with emit-ters polarized transverse to the array (DD-T). (b) Dispersionrelation ω k and decay rate γ k for DD-T with k = 0 . π/d .The interval of [ − k , k ] is shaded. (c) Scaling in N of thedecay rates of the most subradiant states, for DD-T with k dπ = 0 . , . , .
55. The dashed lines are guides of eyes.(d) Amplitudes profile of the first (red) and second (grey)most subradiant states. of DD-T with k = 0 . π/d (dashed)and k = 0 . π/d (solid), N = 100. them with wave numbers. In this sense, people foundthe most subradiant states of DD-T near the extremum k ex = π/d where γ k = 0 [4, 6].In Fig. 1(c), we plot the decay rate of the most sub-radiant state as a function of N for DD-T with threedifferent k . The curve of k = 0 . π/d oscillates but de-scends roughly like N − γ , where γ is the spontaneousemission rate. The curve of k = 0 . π/d shows a perfect N − descending. However, the curve of k = 0 . π/d shows an N − scaling, which is not seen before. Am-plitudes profiles of the first two cases are found to besimilar, but visibly different from those of the third, seeFig. 1(d) for a comparison. As a starter for the rela-tion between α and s , we show in [31] that for DD-Twith k = k (4) ≈ . π/d , k = π/d is an extremum of s = 4, while s = 2 otherwise.To finish the preparations, suppose H eff = H Reeff − iH Imeff where H Reeff and − iH Imeff are the coherent and dissipationparts, respectively. It is simple to prove [31] that thediscrete translation symmetry leads to H Imeff = 12 N Z k ∈ Γ dk π γ k | k i h k | . (3)Therein, “Γ” denotes the the closed subset where γ k = 0.For DD-T, we have Γ = [ − k , k ]. Equation (3) also pro-vides a model of general H Imeff . For our purpose, knowl-edge of γ k is unnecessary. But we require that k ex / ∈ Γ. Eigenstates.
We deliberately digress from the themeand consider the eigenstates of a toy Hamiltonian H R = h I + R X r =1 N − r X i =1 ( h r a † i a i + r + h r a † i + r a i ) (4)where a i ( a † i ) is the annihilation (creation) operator ofthe i th site, and R is the farthest hopping distance. Theleftmost R sites ∂ l = { , , · · · R } and (similarly) therightmost R sites ∂ r constitute the boundary set ∂ l ∪ ∂ r ,of which the projector is denoted by P ∂ . Projector of theremaining sites (bulk) is denoted by P . Then the eigen-state satisfying ( H R − E ) | ψ i = 0 must fulfill the bulkequation P ( H R − E ) | ψ i = 0. If we find out all the so-lutions to the bulk equation of eigenvalue E , which forma linear space, the desired eigenstate must be contained.And it must be the one satisfying also the boundary equa-tion P ∂ ( H R − E ) | ψ i = 0.For the bulk equation, we note that the Bloch state | k i , and its generalizations with k being complex num-bers, satisfy P H R | k i = ω R ( k ) P | k i , thus also the bulkequation if ω R ( k ) = E . Actually, in generic situations allsolutions to the bulk equation (a linear space) is spannedby all such {| k i} . This is mathematically proved by ageneralized Bloch’s theorem [32–34]. To proceed, we in-troduce z = e ikd and re-denote the Bloch state of Eq. (2)by | z i . Then ω R ( k ) becomes H R ( z ) = h + R X r =1 ( h r z r + h r z − r ) . (5)From it we know that H R ( z ) = E , a polynomial equation,has totally 2 R solutions. Suppose k ex is an extremum of ω R ( k ) of order s ( s ≤ R ). Then near z ex = e ik ex , wehave the expansion H R ( z ) = H R ( z ex ) + a s iz ex ) s ( z − z ex ) s + · · · . (6)To choose a reasonable E , note that there are totally N eigenstates in the Brillouin Zone of length 2 π/d , hencethe separation between the neighboring two is O ( N − ).Then states near k ex have eigenvalues E = H R ( z ex ) + a s δ s with δ ∼ N − . Equation (6) implies s solutions to H R ( z ) = E distributed close to z ex : z j ≈ z ex (1 + iδe i πs j ) , j = 1 , · · · s. (7)All solutions are here if s = 2 R (always achievable byadjusting { h r , h r } r of H R , we shall call this special toymodel Hamiltonian H s/ ).For the boundary equation, we need to find coefficients { c j } so that | ψ i = P j c j | z j i satisfies P ∂ ( H R − E ) | ψ i =0. Equivalently, we consider the boundary matrix M de-fined with elements M b,j = h b | ( H R − E ) | z j i , and require { c j } to satisfy P j M b,j c j = 0 for any b ∈ ∂ l ∪ ∂ r . The so-lution exists iff det M = 0. By introducing (cid:15) j , η j ∼ N − so that z ex /z j = 1 + (cid:15) j = (1 + η j ) − , we prove in [31]that the boundary equation is equivalent to s X j =1 c j (cid:15) rj = 0 , s X j =1 c j z N +1 j η rj = 0 , (8)where r = 0 , , · · · , s/ − Scalings.
Now we return to the decay rates and con-sider H s/ − iH Im eff (remember that k ex / ∈ Γ). For eigen-states near k ex , we view H Im eff as a perturbation to H s/ .Thus to the leading order, the decay rate is obtained by γ/ h ψ | H Im eff | ψ i . The validity of this treatment will beseen later. Firstly, for k ∈ Γ, h k | ψ i equals h k | ψ i = O ( 1 N ) s X j =1 c j z j e − ikd − ( z j e − ikd ) N +1 − z j e − ikd . (9)We split the fraction and evaluate separately: X sj =1 c j z j e − ikd − z j e − ikd = X sj =1 c j z − j e ikd −
1= 1 z − ex e ikd − ∞ X n =0 P sj =1 c j (cid:15) nj ( z ex e − ikd − n , (10a) X sj =1 c j ( z j e − ikd ) N +1 − z j e − ikd = e − i ( N +1) kd − z ex e − ikd ∞ X n =0 P sj =1 c j z N +1 j η nj ( e ikd − z ex ) n . (10b)By substituting Eq. (8) into Eq. (10), we see that in theseries of n , the leading contribution can only come from n = R = s/
2. Thus h k | ψ i ∼ N − s/ − and h ψ | H Imeff | ψ i ∼ N × | N − s − | ∼ N − s − . (11)Note that the gaps between eigenstates of H s/ near k ex are O ( N − s ), which is much larger than h ψ | H Imeff | ψ i . Thusthe perturbation treatment is valid and we obtain therelation α = s + 1 for H s/ .However, coherent part H Reeff of a real system such asDD-T does not have finite R , or always has 2 R (cid:29) s .Nevertheless, H s/ well describes the neighbor of k ex ,while we can patch H s/ with a new term H Ω for k ∈ Ω,a closed subset of the Brillouin Zone ( k ex / ∈ Ω). Then H s/ + H Ω can match H Reeff well. H Ω surely has the formsame with Eq. (3), because the latter requires only thediscrete translation symmetry [31]. Then the same calcu-lation shows that H Ω is also a perturbation to the eigen-states of H s/ near k ex . Therefore, we conclude that α = s + 1 is indeed a universal property, as what is al-ready shown in Fig. 1(c).Meanwhile, the eigenstate | ψ i of H s/ near k ex will beperturbed by H Ω and − iH Im eff . Because h ψ ξ | H Im eff | ψ ξ i ∼ N − s − (and so does H Ω ), the primary perturbation for-malism tells us that | ψ i → | ψ i ∝ | ψ i + O ( N − ) | ψ ⊥ i (12)where h ψ | ψ ⊥ i = 0. It leads to the scaling of the infideli-ties that 1 − | h ψ | ψ i | ∼ N − .To verify the infidelities and thus verify the perturba-tion approach, we consider a toy model H R =2 with pos-itive couplings h = h and h = h , and the dispersionrelation ω ( k ) = 2 h cos( kd ) + 2 h cos(2 kd ). k ex = π/d isan extremum of s = 4 if h /h = 4, and a non-degenerate extremum of s = 2 if h /h >
4. Eigenstates of H R =2 near k ex = π/d are analytically derived in [31]. In Fig. 2,we plot the infidelities between the most subradiant stateof DD-T ( k d/π = 0 . , .
55) with the eigenstate of H R =2 ( h /h = 6) that is closest to k ex , and similarly betweenDD-T( k = k (4) ) and H s/ ( h /h = 4). The N − scal-ing is perfectly confirmed for k d/π = k (4) and 0.55.For curves of k d/π = 0 .
3, oscillations are found inboth Figs. 1(d) and 2. This is caused by degeneraciesof k ex = π/d . To see it, we plot the second and forthTaylor’s series coefficients ( a , ) of ω k at k = π/d in theupper panel of Fig. 2. Therein, we find a > a < k < k (4) , hence ω k has three extrema of s = 2, as il-lustrated in the insertion panel of Fig. 2. So k ex = π/d is degenerate to two Bloch states. Their contributionsto | ψ ξ i are enhanced if their wavelengths match certainresonance with the length of the chain ( N d ). We high-light that similar oscillations were also found in subradi-ant bound states [35] and firstly explained in [29]. It is ω k −7 γ Figure 2. Upper panel: The 2nd and 4th Taylor’s expansioncoefficients ( α , ) of DD-T dispersion relation ω k = π/d , as afunction of k . Region of k < k (4) is shaded. Zoom in panel:fine structure of ω k around k = π/d for DD-T with k =0 . π/d . Lower panel: infidelities (log scale) between themost subradiant states of DD-T with k d/π = 0 . , .
55 and k (4) , and the “fundamental mode” eigenstate (near k = π/d )of H with h /h = 6 and h /h = 4, respectively. Thedashed lines are guides of eyes. also seen (but more clearly) in the toy model H R =2 with h /h <
4, which is solved analytically in [31].Additionally, for fixed N , as k approaches to k (4) frombelow, the separation between the three extrema will be O ( N − ). They might provide totally six z j ≈ z ex thuslead to an N − scaling, given that we keep adjusting k according to the increasing N [9]. Radiation surge in topological transition.
In systemssupporting multiple bands, the topological transition isusually companied by band gap closing at some point ofthe Brillouin Zone. If that point is an extremum ( k ex / ∈ Γ) before the gap closes, the transition may decreasesits s -index discontinuously. Although the above theorydirectly covers only the simple arrays with one band, itinspires us that the shift of s would induce a surge ofradiation, if the system is prepared in subradiant stateswith wave number close to k ex .We numerically demonstrate this effect in a dimerizedDD-T depicted in Fig. 3(a), where a unit cell has twoemitters separated by d and the period of the array is d . The dimerized chain has two bands, which are topo-logically trivial (absence of topology protected boundarystates) if d < d where d = d − d , and nontrivial if d > d . An SSH type topological transition occurs atthe critical point d = d [36, 37]. In Fig. 3(b), partof the dispersion relations (for k / ∈ Γ) of the dimerizedDD-T with k = 0 . π/d and d /d = 0 . , . , .
53 areplotted, respectively. We find that k ex = π/d is an ex-tremum of s = 2 if d /d = 0 .
5. But at the critical point d /d = 0 .
5, the band gap closes and the dispersion rela-tion is linearized. Then we focus on the subradiant statesassociated with wave numbers closest to k = π/d , in boththe upper and the lower band. Their decay rates are plot-ted in Fig. 3(c). The N − -scaling is seen for both bandswhen d /d = 0 .
47 (topological trivial) and 0 .
53 (topolog-ically nontrivial), but N − -scaling instead, at the criticalpoint d /d = 0 .
5. So the numerical results confirm theexpected surge of radiation.Next, we consider dimerized chains coupling to an ideal1D waveguide where the Hamiltonian is [38] H = − i Γ X i,j e ik | x i − x j | σ † i σ j . (13)This model is dissipative for only k = ± k . Interest-ingly, being directly inspired by Ref. [29], we notice thatin the subspace of single-excitation states the inverse ofEq. (13), H − , is the original SSH model where the stag-gered nearest-neighbor-couplings are J = 1 / sin( k d )and J = 1 / sin( k d ). The only difference is that thetwo emitters sitting at the two chain ends are dissipa-tive [31]. The critical points of this model are d = d and d = d ± π/k . To be different from DD-T, we con-sider the latter case where the band gap closes at k = 0,an extremum of s = 2 if the system is away from thecritical point. Subradiant states associated with wave ( a ) xd d ( b ) ( c ) N −3 N −1 Figure 3. (a) Dimerized DD-T. (b) The dispersion relationof DD-T with k = 0 . π/d and d /d = 0 . , . , .
53, re-spectively. (c) Decay rates of the subradiant states that areassociated with the wave number closest to k ex = π/d . Thedotted (solid) lines refer to subradiant states of the upper(lower) band. The scaling is changed from N − to N − atthe critical point, N is the number of unit cells. The dashedlines are guides of eyes. numbers close to k ex = 0 have α = 3. However, at thetopological transition, their decay rates become (assum-ing sin( k d ) > γ = Γ N cot( k d ) ln( 1 + sin k d − sin k d ) , (14)for states in both bands [31]. Same with DD-T, this isalso a jump of α from 3 to 1. Conclusions and Discussions.
In this Letter, we haveestablished a universal equality α = s + 1, where s is theorder of an extremum ( k ex ) of the dispersion relation,and α is the exponent of the decay rates scaling ( N − α )of subradiant states with wave numbers close to k ex . Itanswers the question that what determines α , and impliesa surge of radiation of the subradiant states induced byband gap closing. We demonstrated this effect by a jumpfrom the N − scaling to N − that numerically found indimerized emitter arrays coupled to 3D free space modes,and analytically in the case of ideal 1D waveguide.To substantiate the universality, we developed an ap-proach independent to specific Hamiltonians. We believethat our method can be extended to arrays of more com-plex lattices and higher dimensions. In 2D and 3D, geom-etry of the lattices and energy bands will have much morevarieties. Thus a more sophisticated relation might beexpected. For other possible future studies, we note thatthis work reveals the intrinsic connection between subra-diant states and flatter bands (larger s ), as what high-lighted in [29]. Thus applications relying on subradiantstates may refer to systems having flatter bands. Pho-tonic flat bands [39], if being coupled to quantum emit-ters, may provide new controllable platforms for subradi-ant states, and for strongly-correlated manybody physicssuch as the fractional Hall effect that relevant to flatbands [40].The author acknowledges financial support of Euro-pean Unions Horizon 2020 Research and Innovation Pro-gram under Grant Agreement No. 820445 (Quantum In-ternet Alliance). ∗ [email protected][1] P. Weiss, M. O. Ara´ujo, R. Kaiser, and W. Guerin, NewJ. Phys. , 063024 (2018).[2] R. H. Dicke, Phys. Rev. , 99 (1954).[3] H. R. Haakh, S. Faez, and V. Sandoghdar, Phys. Rev.A , 053840 (2016).[4] A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht,H. J. Kimble, and D. E. Chang, Phys. Rev. X , 031024(2017).[5] A. Albrecht, L. Henriet, A. Asenjo-Garcia, P. B. Dieterle,O. Painter, and D. E. Chang, New J. Phys. , 025003(2019).[6] Y.-X. Zhang and K. Mølmer, Phys. Rev. Lett. ,203605 (2019).[7] T. Yu, Y.-X. Zhang, S. Sharma, X. Zhang, Y. M. Blanter,and G. E. W. Bauer, Phys. Rev. Lett. , 107202(2020).[8] F. Dinc, L. E. Hayward, and A. M. Bra´nczyk, “Multi-dimensional super- and subradiance in waveguide quan-tum electrodynamics,” (2020), arXiv:2003.04906 [quant-ph].[9] D. F. Kornovan, N. V. Corzo, J. Laurat, and A. S.Sheremet, Phys. Rev. A , 063832 (2019).[10] A. Asenjo-Garcia, H. J. Kimble, and D. E. Chang, Pro-ceedings of the National Academy of Sciences , 25503(2019).[11] S. D. Jenkins, J. Ruostekoski, N. Papasimakis, S. Savo,and N. I. Zheludev, Phys. Rev. Lett. , 053901 (2017).[12] M. Mirhosseini, E. Kim, X. Zhang, A. Sipahigil, P. B.Dieterle, A. J. Keller, A. Asenjo-Garcia, D. E. Chang,and O. Painter, Nature , 692 (2019).[13] P.-O. Guimond, A. Grankin, D. V. Vasilyev, B. Vermer-sch, and P. Zoller, Phys. Rev. Lett. , 093601 (2019).[14] E. Shahmoon, D. S. Wild, M. D. Lukin, and S. F. Yelin,Phys. Rev. Lett. , 113601 (2017).[15] J. Rui, D. Wei, A. Rubio-Abadal, S. Hollerith, J. Zeiher,D. M. Stamper-Kurn, C. Gross, and I. Bloch, “A subra-diant optical mirror formed by a single structured atomiclayer,” (2020), arXiv:2001.00795 [quant-ph].[16] R. Bekenstein, I. Pikovski, H. Pichler, E. Shahmoon, S. F.Yelin, and M. D. Lukin, Nature Physics , 676 (2020). [17] Y.-X. Zhang, Y. Zhang, and K. Mølmer, Phys. Rev. A , 033821 (2018).[18] Y. Ke, A. V. Poshakinskiy, C. Lee, Y. S. Kivshar, andA. N. Poddubny, Phys. Rev. Lett. , 253601 (2019).[19] N. J. Schilder, C. Sauvan, Y. R. P. Sortais, A. Browaeys,and J.-J. Greffet, Phys. Rev. Lett. , 073403 (2020).[20] R. J. Bettles, M. D. Lee, S. A. Gardiner, andJ. Ruostekoski, “Quantum and nonlinear effects in lighttransmitted through planar atomic arrays,” (2019),arXiv:1907.07030 [quant-ph].[21] J. Cremer, D. Plankensteiner, M. Moreno-Cardoner,L. Ostermann, and H. Ritsch, “Polarization control ofradiation and energy flow in dipole-coupled nanorings,”(2020), arXiv:2004.09861 [quant-ph].[22] G. Facchinetti, S. D. Jenkins, and J. Ruostekoski, Phys.Rev. Lett. , 243601 (2016).[23] M. T. Manzoni, M. Moreno-Cardoner, A. Asenjo-Garcia,J. V. Porto, A. V. Gorshkov, and D. E. Chang, NewJournal of Physics , 083048 (2018).[24] M. Moreno-Cardoner, D. Plankensteiner, L. Ostermann,D. E. Chang, and H. Ritsch, Phys. Rev. A , 023806(2019).[25] J. A. Needham, I. Lesanovsky, and B. Olmos,arXiv:1905.00508.[26] K. E. Ballantine and J. Ruostekoski, Physical ReviewResearch (2020), 10.1103/physrevresearch.2.023086.[27] J. Perczel, J. Borregaard, D. E. Chang, H. Pichler, S. F.Yelin, P. Zoller, and M. D. Lukin, Phys. Rev. Lett. ,023603 (2017).[28] R. J. Bettles, J. c. v. Min´aˇr, C. S. Adams, I. Lesanovsky,and B. Olmos, Phys. Rev. A , 041603 (2017).[29] A. N. Poddubny, Phys. Rev. A , 043845 (2020).[30] H. T. Dung, L. Kn¨oll, and D.-G. Welsch, Phys. Rev. A , 063810 (2002).[31] Supplemental Material .[32] A. Alase, E. Cobanera, G. Ortiz, and L. Viola, Phys.Rev. Lett. , 076804 (2016).[33] E. Cobanera, A. Alase, G. Ortiz, and L. Viola, Journalof Physics A: Mathematical and Theoretical , 195204(2017).[34] A. Alase, E. Cobanera, G. Ortiz, and L. Viola, Phys.Rev. B , 195133 (2017).[35] Y.-X. Zhang, C. Yu, and K. Mølmer, Phys. Rev. Re-search , 013173 (2020).[36] B. X. Wang and C. Y. Zhao, Phys. Rev. A , 023808(2018).[37] S. R. Pocock, X. Xiao, P. A. Huidobro, and V. Giannini,ACS Photonics , 2271 (2018).[38] D. E. Chang, L. Jiang, A. V. Gorshkov, and H. J. Kim-ble, New J. Phys. , 063003 (2012).[39] D. Leykam and S. Flach, APL Photonics , 070901(2018), https://doi.org/10.1063/1.5034365.[40] S. A. Parameswaran, R. Roy, and S. L. Sondhi, ComptesRendus Physique , 816 (2013). upplemental Material to “Universal Scaling of Decay Rates of Emitter Arrays in theSubradiant States” Yu-Xiang Zhang
Niels Bohr Institutet, Københavns Universitet, 2100 København Ø, Denmark (Dated: June 12, 2020)
This Supplemental Material is arranged as following:In Sec. S-I, we show that DD-T has an extremum of s = 4when k = k (4) ; In Sec. S-II, we show that Hamiltonianshaving the discrete translation symmetry, including H Im eff and H Ω , can be written into the form of Eq. (3) of themain text; In Sec. S-III, we derive the equivalent forms ofthe boundary equation and Eq. (8) of the main text; InSec. S-IV, we analytically derive the eigenstates of the toymodel H R =2 for cases of h = 4 h ( s = 4) and h > h ( s = 2, k ex is non-degenerate), and h < h ( s = 2, k ex is degenerate). In Sec. S-V, we show the inverse of H and analytically derive Eq. (13) of the main text. S-I. DD-T: THE s -INDEX OF k ex = π/d The dispersion relation of DD-T is given by [S1] ω k = 3 γ X ε = ± X ξ =1 i ( ik d ) ξ Li ξ [ e i ( k + εk ) d ] , where Li ξ denotes the polylogarithm of order ξ Li ξ ( z ) = P ∞ n =1 z n n − ξ . The series does not converge if ξ <
1. Inthis case, we can use Li ( z ) = − Ln(1 − z ) and calculatethem via the relation ddz Li ξ ( z ) = 1 z Li ξ − ( z ) . For DD-T, the parity symmetry ( ω k = ω − k ) makes allthe odd order derivatives of ω k at k = π/d vanish. Thenthe second order derivative is found to be ∂ kd ω kd = π = 32 k d (cid:20) ln(2 cos k d k d k d − ( k d k d ) (cid:21) . Numerical calculation shows that it vanishes at k = k (4) ≈ . π/d . S-II. DERIVATION OF EQ. (3)
Consider H = P Ni,j =1 H ij | i i h j | , where H ij dependsonly on i − j (discrete translation symmetry). The dis-persion relation ω k is obtained from the Fourier transfor-mation of H ij H ij = Z dk/ (2 π ) ω k e ik ( x i − x j ) (S1) so that H = N X i,j =1 Z B.Z. dk π ω k e ik ( x i − x j ) | i i h j | = N Z B.Z. dk π ω k | k i h k | , (S2)where “B.Z.” refers to the first Brillouin zone. In themodel of H Imeff , the corresponding “dispersion relation” isjust γ k / k ∈ Γ, a subset of the Brillouin zone.For DD-T, we showed in Ref. [S2] that the effectiveHamiltonian of DD-T can be written as H eff = − i γ k Z k d ˜ k π ρ + (˜ k ) N X m,n =1 e i ˜ k | z m − z n | σ † m σ n − γ k Z + ∞ d ˜ k π ρ − (˜ k ) N X m,n =1 e − ˜ k | z m − z n | σ † m σ n . where ρ ± (˜ k ) = π (1 ± ˜ k /k ). Therefore in the single-excitation sector, we have H Imeff = N γ k Z k − k d ˜ k π ρ + (˜ k ) | ˜ k i h ˜ k | . So DD-T Hamiltonian can be modeled by Eq. (3) of themain text.
S-III. EQUIVALENT FORMS OF THEBOUNDARY MATRIX
The boundary equation P ∂ ( H R − E ) | z j i = 0 can beformulated through the boundary matrix M . For conve-nience, here we use the notation that | z j i = N X m =1 z mj | m i . Then the element M b,j ≡ h b | ( H R − E ) | z j i reads M b,j = ( h − E ) z bj + R X r =1 ( θ b + r h r z b + rj + θ b − r h r z b − rj ) . (S3)Therein, we introduced the notation θ x , which equals 1if x ∈ ∂ l ∪ ∂ r and equals 0 otherwise. The bulk equation H ( z j ) = E is written as h + R X r =1 ( h r z rj + h r z − rj ) = E. a r X i v : . [ qu a n t - ph ] J un By substituting the above equation into Eq. (S3), we caneliminate h − E and obtain M b,j = z bj R X r =1 [(1 − θ b + r ) h r z rj + (1 − θ b − r ) h r z − rj ]To write them explicitly, for b ∈ ∂ l = { , , · · · R } , wehave M R,j = h R ,M R − ,j = h R z − j + h R − , · · · M ,j = h R z − ( R − j + · · · + h z − j + h . (S4)Expressions for b ∈ ∂ r = { N, N − , · · · , N − R + 1 } willbe given below.We need to find the kernel of the boundary matrix M , i.e., to find { c j } j satisfying P j M b,j c j = 0. For thispurpose, we are free to manipulate M by row operations:to divide one row by some number, and to sum one rowwith another row multiplied by some factor. The firststep is M R,j → M R,j /h R = 1 . Next, we work on the second row: M R − ,j → ( M R − ,j − h R − M R,j ) /h R = z − j . Similarly for the remaining rows. Finally, we obtain anormal form of M b,j ( b ∈ ∂ l ): M R,j = 1 ,M R − ,j = z − j , · · · M ,j = z − ( R − j . (S5)For M b,j ( b ∈ ∂ r ), we have the normal form M N − R +1 ,j = z N +1 j ,M N − R +2 ,j = z N +2 j , · · · M N,j = z N + Rj . (S6)Based on that, it is more convenient to denote rows of b ∈ ∂ l as M L and those of b ∈ ∂ r as M R , and relabel thematrix elements as M Lr,j = z − rj , M Rr,j = z N + rj , (S7)where r = 1 , , · · · R . The equation P j M b,j c j = 0 forevery b is equivalent to require X j z − rj c j = 0 , X j z N + rj c j = 0 . (S8) for every r ∈ { , , · · · R } .The boundary matrix elements can be reformulated by (cid:15) j and η j . Firstly we multiply the r-th row of M L by z r − so that M Lr,j = (1 + (cid:15) j ) r − , or M Lr +1 ,j = r X k =0 C kr (cid:15) kj , where C kr is the binomial coefficient. Then we sequen-tially apply the row operations M Lr +1 ,j → M Lr +1 ,j − r − X m =0 C mr M Lm,j by the order of r = 1 , , · · · , to reshape the rows into M L = · · · (cid:15) · · · (cid:15) s ... ... ... (cid:15) R − · · · (cid:15) R − s R × s (S9)As to M R , firstly we multiply M Rr +1 ,j by z − r ex so that M Rr +1 ,j = z N +1 j (1 + η j ) r = z N +1 j r X k =0 C kr η kj . Following the similar row manipulations, we have M R = z N +11 · · · z N +1 s z N +11 η · · · z N +1 s η s ... ... ... z N +11 η R − · · · z N +1 s η R − s R × s (S10)Then we obtain Eq. (8) of the main text. S-IV. ANALYTICAL RESULTS OF THE TOYMODEL H R =2 We consider the case of s = 2 and s = 4 in the toymodel ( R = 2) of the Hamiltonian H = h N − X i =1 a † i a i +1 + h N − X i =1 a † i a i +2 . (S11)Therein, we assume h and h are both positive numbers.The dispersion relation of this model is given by ω ( k ) = 2 h cos( k ) + 2 h cos(2 k ) . The point k = π/d is an extremum of s = 4 if h =4 h (2 R = s , essential solution), and a non-degenerateextremum of s = 2 if h > h .To derive the eigenstates near k ex , firstly we study thebulk equation H ( z ) = h ( z + 1 z ) + h ( z + 1 z ) = E. (S12)Therein, we suppose E = ω ( k ex ) + e with e > e (cid:28) h . The equation H ( z ) = E leads to z + 1 z = − h h ± h p ( h − h ) + 4 h e. (S13)Suppose we write the boundary matrix M in the normalform of Eqs. (S5) and (S6). Then the determinant of M is given bydet M ∝ X j 5. Then we have − c c = e − ζπ + ( − ξ e ζπ − ( − ξ . (S17)If we let c = e ζπ + ( − ξ , the other two coefficients willbe given by c = ( − ξ − sinh( ζπ ) + i cosh( ζπ ) ,c = ( − ξ − sinh( ζπ ) − i cosh( ζπ ) . (S18)Note that the four coefficients have comparable magni-tudes. S-IV.B. Case of 2R > s = > ) We expand the square root in Eq. (S13) to first orderof e and get p ( h − h ) + 4 h e = | h − h | + 2 h e | h − h | . We use the shorthand expression h /h − a ( a > e/ ( ah ) = δ . Then we have z + 1 /z = − δ, ˜ z + 1 / ˜ z = − ( a + 2) − δ. The four solutions are z , and ˜ z , [these two solutionshave magnitudes different with 1 by O (1), the corre-sponding | ˜ z , i are boundary states]: z = 1 /z ≈ − i √ δ ≈ − e − i √ δ , ˜ z = 1 / ˜ z − ≈ − xe βδ , where β = a +2 + √ a +4 a and x = a + 1. Substitutingthem into Eq. (S14) leads to √ δ ( N + 2 − x + 1 x − ≈ ξπ (S19)where ξ is integer. The above approximation is exponen-tially precise [up to O ( x − N )]. Note thatlim x (cid:29) √ δ = ξπN + 1 , (S20)is the essential solution of s = 2 (obtained in the toymodel of R = 1). Solutions given by Eq. (S19) devi-ate from the essential solution by O ( N − ). The relativedeviations of (cid:15) j and η j are thus O ( N − ).For the superposition coefficients, we have˜ c ˜ c = ( − N + ξ ˜ z N +13 , (S21)which shows the balance between them [note that thenorm of | ˜ z i scales like O ( | ˜ z | N +1 ) while that of | ˜ z i scales as O (1)]. Coefficients { c j } j are c = − i − x − i √ δ + βδ √ δ ˜ c ,c = i − x + i √ δ + βδ √ δ ˜ c . (S22)where we have neglect the exponentially small contri-bution from ˜ c . So we see that the eigenstates satisfyEq. (12) of the main text. S-IV.C. Case of k ex is degenerate ( h < h ) Now k ex = π/d is a maximum. Thus we consider eigen-values E = ω ( k ex ) − e where e > 0. Similar to previoussections, we define 4 − h /h = a and e/ ( ah ) = δ so that z + 1 /z = − δ, z + 1 /z = − a − δ. The solutions are z = 1 /z ≈ − e i √ δ , and˜ z = 1˜ z = − a − δ i p a − a + ( a − δ √ a − a ) . Compared with the previous section, here | ˜ z | = | ˜ z | = 1.Thus we can introduce θ so that ˜ z = e iθ . The depen-dence to δ is expressed as θ = θ + δ , where δ ∼ δ . Thendet M = 0 leads to √ δ = ξπN + 2 (1 + η ) , where ξ = 1 , , · · · and the higher order correction η = − tan( θ / N + 2 1 − ( − ξ + N cos[ θ ( N + 2)]( − ξ + N sin[ θ ( N + 2)] . (S23)Comparing the above equation with Eq. (S19). We seethat now η has trigonometric functions of θ ( N + 2),which lead to the oscillations depending on N .As to the superposition coefficients, suppose c = 1.Then we have c = − e iφ ≈ − − iφ , where φ ≈ ηξπ ( − N + ξ sin θ + √ δ sin[ θ ( N + 2)]( − N + ξ sin θ + sin[ θ ( N + 2)] + sin[ θ ( N + 1)] . Coefficients of | ˜ z i are˜ c ≈ √ δ sin θ − e − iθ θ φ and ˜ c ≈ − ˜ c + iφ . Thus the relation c = O ( N )˜ c is still valid and the eigenstates satisfy Eq. (12) of themain text. S-V. DIMERIZED IDEAL 1D WAVEGUIDEMODEL We consider a chain of N unit cells, hence totally 2 N emitters. The inverse of the 2 N × N effective Hamilto-nian is expressed as H − = 1Γ a J · · · J e J · · · J e J · · · J e J · · ·· · · · · · (S24)where the parameters are a = i − cot( k d ) , J = 1sin( k d ) ,J = 1sin( k d ) , e = − sin( k d )sin( k d ) sin( k d ) . Equation (S24) describes the SSH model where the stag-gered couplings between nearest neighbors are J and J ,and the on-site energy is e . The two chain ends haveshifted on-site energy and also dissipation.In an infinite chain, the eigenstates are Bloch states inthe form of | k, u k i = 1 √ N N X m =1 e ikz m u k · ~σ † m | g a g b · · · g Na g Nb i , where ~σ m = ( σ m,a , σ m,b ), indices “a(b)” refer to the left(right) emitter of a unit cell, and u k is a normalized vec-tor that determines the amplitudes and relative phase ofthe excitations intra a unit cell.The dispersion relation is expressed as ω ± ( k ) = Γ / kd ) − cos( k d ) (cid:20) sin( k d ) ± p J − + ( J ) − + 2( JJ ) − cos( kd ) (cid:21) . (S25)Note that the band ω + diverges at k = ± k (we assumesin( k d ) > ω − ( ± k ) = − Γ k d ) sin( k d )sin( k d ) . (S26)The intra-cell “spin” of the two bands are | u + k i = 1 √ (cid:18) e iφ k (cid:19) , | u − k i = 1 √ (cid:18) − e iφ k (cid:19) . where φ k is given bytan φ k = sin( kd )cos( kd ) + J/J . We see from Eq. (S25) that the band gap closes at k = 0if sin( k d ) = − sin( k d ), and at k = π/d if sin( k d ) =sin( k d ). In the following, we focus only on the formercase. The latter case can be studied in the same way.At the critical point sin( k d ) = − sin( k d ), the dis-persion relation is given by ω ( k ) = Γ / kd ) − cos( k d ) × (cid:20) sin( k d ) ± k d ) sin( kd (cid:21) . (S27)It is linear near k = 0. To compare, k = 0 is an extremumof s = 2 away from the critical point.The intra-cell “spin” shows a discontinuity, e.g., for theupper band we have | u + k> i = 1 √ (cid:18) ie i kd (cid:19) , | u + k< i = − √ (cid:18) − ie i kd (cid:19) . Next, we derive the subradiant states near k = 0 atthe critical point. For abbreviations, we shall use | k i asa short hand of an ansatz | k, n k i , and use | k ± i to denote | k, u ± k i . For finite chains, we have H | k i = ω (cid:15)k h k (cid:15) | k i | k (cid:15) i + g k | k , + i − h k |− k , + i , where (cid:15) = ± and g k = 1 i √ e i ( k − k ) z − e i ( k − k ) d h u + k | n k i , h k = 1 i √ e i ( k + k )( z N + d ) − e i ( k + k ) d h u + − k | n k i . The eigenstates are superpositions of two Bloch states |± k i which satisfy g k h − k = g − k h k . (S29) S-V.A. the upper band Equation (S29) is expanded as e − iNdk sin k + k d sin k − k d = 1 − sin( k d + kd )1 + sin( k d − kd ) . (S30)For those k ≈ 0, we substitute the ansatz k = (cid:15) πN d (1 + 1 N δ )into Eq. (S30) and obtain e − i π(cid:15) = 1 − sin k d k d ,iδ = cot( k d k d ) . The solution is not unique and we have (cid:15) ξ = ξ + i π ln 1 − sin k d k d , where ξ = 1 , , · · · . The imaginary parts of (cid:15) ξ measuresthe biased mixture of the two components of |± k + i .Substituting the solutions of k into ω + ( k ), we obtainthe frequency shift and the decay rates ω + ( ξ ) = Γ / − cos k d (sin k d + ξ πN sin k d ) , (S32a) γ + = − Γ N cot( k d ) ln( 1 − sin k d k d ) . (S32b)Note that γ + is uniform for all ξ , and γ + ∝ N − . S-V.B. the lower band Now Eq. (S29) yields e − iNdk sin k + k d sin k − k d = 1 + sin( k d + kd )1 − sin( k d − kd ) . Therefore the solutions would be (cid:15) ξ = ξ + i π ln 1 + sin k d − sin k d . And the eigenvalues are given by ω − ( ξ ) = Γ / − cos k d (sin k d − ξ πN sin k d ) , (S33a) γ − = Γ N cot( k d ) ln( 1 + sin k d − sin k d ) . (S33b) [S1] A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht,H. J. Kimble, and D. E. Chang, Phys. Rev. X , 031024(2017). [S2] Y.-X. Zhang and K. Mølmer, Phys. Rev. Lett.122