SQUID Metamaterials on a Lieb lattice: From flat-band to nonlinear localization

TL;DR

This paper investigates the linear and nonlinear localization phenomena in SQUID metamaterials arranged on a Lieb lattice, revealing flat-band and discrete breather states, and analyzing their dependence on initial excitation and amplitude.

Contribution

It introduces a detailed analysis of flat-band and nonlinear localized states in SQUID Lieb lattices, combining linear spectrum derivation with numerical simulations of localization regimes.

Findings

Flat-band localized states occur in nearly linear regimes.

Nonlinear discrete breather states form at higher amplitudes.

Localization depends on initial excitation position and amplitude.

Abstract

The dynamic equations for the fluxes through the SQUIDs that form a two-dimensional metamamaterial on a Lieb lattice are derived, and then linearized around zero flux to obtain the linear frequency spectrum according to the standard procedure. That spectrum, due to the Lieb lattice geometry, possesses a frequency band structure exhibiting two characteristic features; two dispersive bands, which form a Dirac cone at the corners of the first Brillouin zone, and a flat band crossing the Dirac points. It is demonstrated numerically that localized states can be excited in the system when it is initialized with single-site excitations; depending on the amplitude of those initial states, the localization is either due to the flat-band or to nonlinear effects. Flat-band localized states are formed in the nearly linear regime, while localized excitations of the discrete breather type are formed in the nonlinear regime. These two regimes are separated by an intermediate turbulent regime for which no localization is observed. Notably, initial single-site excitations of only edge SQUIDs of a unit cell may end-up in flat-band localized states; no such states are formed for initial single-site excitations of a corner SQUID of a unit cell. The degree of localization of the resulting states is in any case quantified using well-established measures such as the energetic participation ratio and the second moment.

Equations

$${\Phi }_{n,m}^A={\Phi }_{ext}+L\{I_{n,m}^A+{\lambda }_x\left[I_{n-1,m}^B+I_{n,m}^B\right]$$

$$+{\lambda }_y\left[I_{n,m-1}^C+I_{n,m}^C\right]\},$$

$${\Phi }_{n,m}^B={\Phi }_{ext}+L\left\{I_{n,m}^B+{\lambda }_x\left[I_{n,m}^A+I_{n+1,m}^A\right]\right\},$$

$${\Phi }_{n,m}^C={\Phi }_{ext}+L\left\{I_{n,m}^C+{\lambda }_y\left[I_{n,m}^A+I_{n,m+1}^A\right]\right\},$$

$$-I_{n,m}^k=C\frac{d^2{\Phi }_{n,m}^k}{d{t}^2}+\frac{1}{R}\frac{d{\Phi }_{n,m}^k}{dt}+I_c\sin \left(2\pi \frac{{\Phi }_{n,m}^k}{{\Phi }_0}\right),$$

$$L{I}_{n,m}^A={\Phi }_{n,m}^A-{\lambda }_x\left({\Phi }_{n,m}^B+{\Phi }_{n-1,m}^B\right)$$

$$-{\lambda }_y\left({\Phi }_{n,m}^C+{\Phi }_{n,m-1}^C\right)-{\Phi }_{eff}^A,$$

$$L{I}_{n,m}^B={\Phi }_{n,m}^B-{\lambda }_x\left({\Phi }_{n,m}^A+{\Phi }_{n+1,m}^A\right)-{\Phi }_{eff}^B,$$

$$L{I}_{n,m}^C={\Phi }_{n,m}^C-{\lambda }_y\left({\Phi }_{n,m}^A+{\Phi }_{n,m+1}^A\right)-{\Phi }_{eff}^C,$$

$$LC\frac{d^2{\Phi }_{n,m}^A}{d{t}^2}+\frac{L}{R}\frac{d{\Phi }_{n,m}^A}{dt}+L{I}_c\sin \left(2\pi \frac{{\Phi }_{n,m}^A}{{\Phi }_0}\right)+{\Phi }_{n,m}^A$$

$$={\lambda }_x\left({\Phi }_{n,m}^B+{\Phi }_{n-1,m}^B\right)+{\lambda }_y\left({\Phi }_{n,m}^C+{\Phi }_{n,m-1}^C\right)+{\Phi }_{eff}^A,$$

$$LC\frac{d^2{\Phi }_{n,m}^B}{d{t}^2}+\frac{L}{R}\frac{d{\Phi }_{n,m}^B}{dt}+L{I}_c\sin \left(2\pi \frac{{\Phi }_{n,m}^B}{{\Phi }_0}\right)+{\Phi }_{n,m}^B$$

$$={\lambda }_x\left({\Phi }_{n,m}^A+{\Phi }_{n+1,m}^A\right)+{\Phi }_{eff}^B,$$

$$LC\frac{d^2{\Phi }_{n,m}^C}{d{t}^2}+\frac{L}{R}\frac{d{\Phi }_{n,m}^C}{dt}+L{I}_c\sin \left(2\pi \frac{{\Phi }_{n,m}^C}{{\Phi }_0}\right)+{\Phi }_{n,m}^C$$

$$={\lambda }_y\left({\Phi }_{n,m}^A+{\Phi }_{n,m+1}^A\right)+{\Phi }_{eff}^C.$$

$$\tau ={\omega }_{LC}t,{\varphi }_{n,m}^k=\frac{{\Phi }_{n,m}^k}{{\Phi }_0},{\varphi }_{ext}=\frac{{\Phi }_{ext}}{{\Phi }_0},$$

$${\ddot{\varphi }}_{n,m}^A+\gamma {\dot{\varphi }}_{n,m}^A+\beta \sin \left(2\pi {\varphi }_{n,m}^A\right)+{\varphi }_{n,m}^A$$

$$={\lambda }_x\left({\varphi }_{n,m}^B+{\varphi }_{n-1,m}^B\right)+{\lambda }_y\left({\varphi }_{n,m}^C+{\varphi }_{n,m-1}^C\right)+{\varphi }_{eff}^A,$$

$${\ddot{\varphi }}_{n,m}^B+\gamma {\dot{\varphi }}_{n,m}^B+\beta \sin \left(2\pi {\varphi }_{n,m}^B\right)+{\varphi }_{n,m}^B$$

$$={\lambda }_x\left({\varphi }_{n,m}^A+{\varphi }_{n+1,m}^A\right)+{\varphi }_{eff}^B,$$

$${\ddot{\varphi }}_{n,m}^C+\gamma {\varphi }_{n,m}^C+\beta \sin \left(2\pi {\varphi }_{n,m}^C\right)+{\varphi }_{n,m}^C$$

$$={\lambda }_y\left({\varphi }_{n,m}^A+{\varphi }_{n,m+1}^A\right)+{\varphi }_{eff}^C,$$

$$\beta =\frac{L{I}_c}{{\Phi }_0}\mathrm{and}\gamma ={\omega }_{LC}\frac{L}{R},$$

$${\ddot{\varphi }}_{n,m}^A+{\Omega }_{SQ}^2{\varphi }_{n,m}^A={\lambda }_x\left({\varphi }_{n,m}^B+{\varphi }_{n-1,m}^B\right)$$

$$+{\lambda }_y\left({\varphi }_{n,m}^C+{\varphi }_{n,m-1}^C\right),$$

$$\ddot{{\varphi }_{n,m}^B}+{\Omega }_{SQ}^2{\varphi }_{n,m}^B={\lambda }_x\left({\varphi }_{n,m}^A+{\varphi }_{n+1,m}^A\right),$$

$${\ddot{\varphi }}_{n,m}^C+{\Omega }_{SQ}^2{\varphi }_{n,m}^C={\lambda }_y\left({\varphi }_{n,m}^A+{\varphi }_{n,m+1}^A\right),$$

$${\varphi }_{n,m}^k={\mathcal{F}}_k\exp [i(\Omega \tau -{\kappa }_{x}n-{\kappa }_{y}m)],$$

$$\left({\Omega }_{SQ}^2-{\Omega }^2\right){\mathcal{F}}_A-{\lambda }_x\left(1+e^{+i{\kappa }_x}\right){\mathcal{F}}_B$$

$$-{\lambda }_y\left(1+e^{+i{\kappa }_y}\right){\mathcal{F}}_C=0,$$

$$-{\lambda }_x\left(1+e^{-i{\kappa }_x}\right){\mathcal{F}}_A+\left({\Omega }_{SQ}^2-{\Omega }^2\right){\mathcal{F}}_B=0,$$

$$-{\lambda }_y\left(1+e^{-i{\kappa }_y}\right){\mathcal{F}}_A+\left({\Omega }_{SQ}^2-{\Omega }^2\right){\mathcal{F}}_C=0.$$

$$\mathcal{D}=\left({\Omega }_{SQ}^2-{\Omega }^2\right)\{{\left({\Omega }_{SQ}^2-{\Omega }^2\right)}^2$$

$$-4\left[{\lambda }_x^2{\cos }^2\left(\frac{{\kappa }_x}{2}\right)+{\lambda }_y^2{\cos }^2\left(\frac{{\kappa }_y}{2}\right)\right]\}=0.$$

$${\Omega }_{\kappa }$$

$${\Omega }_{\stackrel{max}{\min }}=\sqrt{{\Omega }_{SQ}^2\pm 2\sqrt{{\lambda }_x^2+{\lambda }_y^2}}.$$

$$H=\sum _{n,m}{H}_{n,m},$$

$$H_{n,m}=\sum _k\left\{\frac{\pi }{\beta }\left[{\left(q_{n,m}^k\right)}^2+{\left({\varphi }_{n,m}^k\right)}^2\right]-\cos \left(2\pi {\varphi }_{n,m}^k\right)\right\}$$

$$-\frac{\pi }{\beta }\{{\lambda }_x[{\varphi }_{n,m}^A{\varphi }_{n-1,m}^B+2{\varphi }_{n,m}^A{\varphi }_{n,m}^B+{\varphi }_{n,m}^B{\varphi }_{n+1,m}^A]$$

$$+{\lambda }_y[{\varphi }_{n,m}^A{\varphi }_{n,m-1}^C+2{\varphi }_{n,m}^A{\varphi }_{n,m}^C+{\varphi }_{n,m}^C{\varphi }_{n,m+1}^A]\},$$

$${\vec{\varphi }}_{n+\frac{1}{2}}^k={\vec{\varphi }}_n^k+\frac{h}{2}{\vec{q}}_n^k,$$

$${\vec{q}}_{n+1}^k={\vec{q}}_n^k-h{H}_{{\vec{\varphi }}^k}({\vec{\varphi }}_{n+\frac{1}{2}}^k),$$

$${\vec{\varphi }}_{n+1}^k={\vec{\varphi }}_{n+\frac{1}{2}}^k+\frac{h}{2}{\vec{q}}_{n+1}^k,$$

Full text

SQUID metamaterials, Lieb lattice, flat dispersion bands, localized flat-band states, Hamiltonian systems, Discrete breathers

I Introduction

Considerable research effort has focused the last two decades in the investigation and developement of artificial mediums or metamaterials, which exhibit properties not found in natural materials Smith2004 ; Yen2004 ; Linden2004 ; Linden2006 ; Shalaev2007 ; Litchinitser2008 . After the development of active, tunable, and nonlinear metamaterials Boardman2010 ; Lapine2014 , those artificial mediums are expected to have a strong impact across the entire range of technologies where electromagnetic radiation is used. Moreover, they may provide a flexible platform for modeling and mimicking fundamental physical effects Zheludev2010 ; Zheludev2011 ; Zheludev2012 . An imporant class of metamaterials is that of superconducting ones Anlage2011 ; Jung2014 , and in particular those comprising Superconducting QUantum Interference Devices (SQUIDs). The idea of a metamaterial consisting of SQUIDs was theoretically introduced about a decade ago both in the quantum Du2006 and the classical regimes Lazarides2007 .

The simplest version of a SQUID consists of a superconducting ring interrupted by a Josephson junction Josephson1962 , as shown schematically in Fig. 1. The SQUIDs are highly nonlinear devices, exhibiting strong resonant response to applied magnetic fields. SQUID metamaterials in one and two dimensions have been realized and investigated in the laboratory, and they were found to exhibit novel properties such as negative diamagnetic permeability Jung2013 ; Butz2013a , broad-band tunability Butz2013a ; Trepanier2013 , self-induced broad-band transparency Zhang2015 , as well as dynamic multistability and switching Jung2014b , among others. Some of these properties, i.e., dynamic multistability and tunability, have been also revealed in numerical simulations Lazarides2013b ; Tsironis2014b . Moreover, nonlinear localization Lazarides2008a and the emergence of counter-intuitive dynamic states referred to as chimera states in current literature Lazarides2015b ; Hizanidis2016a have been demonstrated numerically in SQUID metamaterial models.

The notion of metamaterials implies the freedom to engineer not only the properties of the individual ”particles” or devices which play the role of ”atoms” in an artificial medium, but also their arrangement in space, i.e., the type of the lattice. Remarkably, some specific lattice geometries such as those of Lieb or Kagomé lattices give rise to novel and potentially useful band structures. The former is a square-depleted (line-centered tetragonal) lattice, described by three sites in a square unit cell as illustrated in Fig. 2. It is characterized by a band structure featuring Dirac cones intersected by a topological flat band. Localization on flat-bands has been extensively investigated in relatively simple lattice models Danieli2015 ; Khomeriki2016 , even in the presence of disorder Leykam2013 ; Leykam2017 . Superpositions of flat-band modes and their stability have been also investigated in rhombic nonlinear optical waveguide arrays Maimistov2017 . The Lieb lattice was first introduced in the context of photonics in Ref. Leykam2012 . Recently, photonic Lieb lattices have been experimentally realized and the existence of localized flat-band modes has been reported Mukherjee2015a ; Vicencio2015 . The world of electronic flat-band systems has been reviewed in a recent article LiuZheng2014 . Moreover, electronic Lieb lattices have been experimentally realized and characterized Slot2017 . Here, a SQUID metamaterial on a Lieb lattice is considered, in which each site is occupied by a SQUID. In each unit cell, two of the SQUIDs (indicated in red and blue) are neighbored by two other SQUIDs. The third SQUID in the unit cell (black) has four neighbors. In what follows, these SQUIDs will be referred to as edge SQUIDs (red and blue) and corner (black) SQUID, respectively.

In the following, the dynamical equations for the fluxes through the SQUIDs and the linear frequency spectrum are obtained for a SQUID Lieb metamaterial (SLiMM). Using numerical simulations, the generation of localized flat band states when a single edge SQUID is initially excited at low amplitude, is demonstrated. No flat-band localization is observed when single corner SQUIDs are initially excited at low amplitudes, in agreement with the experiments in optical Lieb lattices Guzman-Silva2014 ; Mukherjee2015a ; Vicencio2015 . For high-amplitude initial excitations of either a corner or an edge SQUID, nonlinear localization of the discrete breather type is observed Flach2008a . The cross-over between flat-band and nonlinear localization is explored; the two regimes are clearly separated by an intermediate, no-localization regime. Thus, flat-band localized states cannot be continuated into the nonlinearly localized ones of the discrete breather or discrete soliton type, as it has been demonstrated for discrete nonlinear Schrödinger type models of various flat-band lattices and ribbons Vicencio2013 ; Leykam2013 ; Johansson2015 ; Belicev2015 ; Lopez-Gonzalez2016 ; Gligoric2016 .

II Flux Dynamics

Consider the Lieb lattice of Fig. 2, in which each site is occupied by a SQUID. That SLiMM can be regarded as the combination of three sublattices colored as blue, red, and black. All the SQUIDs are identical, and they are magnetically coupled to their nearest neighbors through their mutual inductances. In order to derive the dynamic equations for the fluxes through the SQUIDs of the SLiMM, we first write the flux-balance relations for all SQUIDs

[TABLE]

where $I_{n,m}^{k}$ is the current in the SQUID of the $(n,m)$th unit cell of kind $k$ ($k=A,B,C$), $\Phi_{ext}$ is the applied (external) flux, and $\lambda_{x}=M_{x}/L$ ($\lambda_{y}=M_{y}/L$) is the coupling coefficient along the horizontal (vertical) direction, with $M_{x}$ ($M_{y}$) being the corresponding mutual inductance between neighboring SQUIDs and $L$ the self-inductance of each SQUID. The current in each SQUID is given by the resistively and capacitively shunted junction (RCSJ) model Likharev1986 , as

[TABLE]

where $R$ is the quasiparticle resistance through the Josephson junction of each SQUID, $C$ is the capacitance of each SQUID, and $I_{c}$ is the critical current of the Josephson junction of each SQUID. Then Eqs. (1) are inverted to obtain the currents $I_{n,m}^{k}$ as functions of the fluxes $\Phi_{n,m}^{k}$. By substitution of the obtained currents back into Eqs. (1), and neglecting all the terms which are proportional to $\lambda_{x}^{a}\lambda_{y}^{b}$ with $a+b>1$, we get

[TABLE]

where $\Phi_{eff}^{A}=[1-2(\lambda_{x}+\lambda_{y})]\Phi_{ext}$, $\Phi_{eff}^{B}=(1-2\lambda_{x})\Phi_{ext}$, and $\Phi_{eff}^{C}=(1-2\lambda_{y})\Phi_{ext}$ are the ”effective” external fluxes. Combining Eqs. (2) and (3) we get

[TABLE]

Using the relations

[TABLE]

where $\omega_{LC}={1}/{\sqrt{LC}}$ is the inductive-capacitive SQUID frequency, the dynamic equations for the fluxes through the SQUIDs can be written in the normalized form

[TABLE]

where

[TABLE]

is the SQUID parameter and the loss coefficient, respectively, and $\phi_{eff}^{k}$ are the normalized effective fluxes. The overdots on $\phi_{n,m}^{k}$ denote differentiation with respect to the normalized temporal variable $\tau$. The values of the fluxes through the SQUIDs $\phi_{n,m}^{k}$ generally depend on $k$. Suppose that $\gamma=0$ and $\phi_{ext}=0$, and that Eqs. (II) are initialized with a low amplitude homogeneous excitation, i.e., with $\phi_{n,m}^{k}=c$ for any $n$, $m$, and $k$ ($c\ll 1$ is a constant). After integrating Eqs. (II) in time assuming periodic boundary conditions, at the steady state, the fluxes through the SQUIDs of the same kind will be the same. However, the fluxes through the SQUIDs of different kind will be different. This is due to the Lieb lattice geometry and the (generally) different values of $\lambda_{x}$ and $\lambda_{y}$, since the flux through a particular SQUID of the SLiMM depends not only on the self-induced one, but also on the fluxes from the SQUIDs to which that particular SQUID is coupled (four for $A$ SQUIDs and two for $B$ and $C$ SQUIDs). Moreover, the coupling between SQUIDs is proportional to the coefficients $\lambda_{x}$ or $\lambda_{y}$. Note that for isotropic coupling, $\lambda_{x}=\lambda_{y}$, the fluxes through the SQUIDs of kind $B$ and $C$ are the same but different than those through the SQUIDs of kind $A$ ($\phi_{n,m}^{B}=\phi_{n,m}^{C}\neq\phi_{n,m}^{A}$).

In the following, we are concerned about energy-conserving SLiMMs, i.e., about the Hamiltonian version of SQUID Lieb metamaterials, and thus we set $\gamma=0$ and $\phi_{ext}=0$ into Eqs. (II).

III Linear Frequency Spectrum

Without losses and driving forces, Eqs. (II) are linearized using the relation $\beta\,\sin\left(2\pi\phi_{n,m}^{k}\right)\simeq\beta_{L}\,\phi_{n,m}^{k}$, where $\beta_{L}=2\pi\beta$. Thus we get

[TABLE]

where $\Omega_{SQ}=\sqrt{1+\beta_{L}}$ is the resonance frequency of individual SQUIDs in the linear limit. In order to obtain the linear frequency spectrum, we substitute into the linearized Eqs. (8) the plane wave solution

[TABLE]

where $\kappa_{x}$ and $\kappa_{y}$ are the $x$ and $y$ components of the two-dimensional, normalized wavevector $\bf\kappa$, and $\Omega=\omega/\omega_{LC}$ is the normalized frequency. After some calculations we get

[TABLE]

In order to obtain nontrivial solutions for the amplitudes ${\cal F}_{k}$ of the stationary problem Eqs. (10), its determinant ${\cal D}$ should be equal to zero, i.e.,

[TABLE]

Solving Eq. (III) for $\Omega\equiv\Omega_{\bf\kappa}$, we get

[TABLE]

where only positive frequencies are considered. Eqs. (12) and (13) provide the linear frequency spectrum of the SLiMM. Thus, the Lieb lattice geometry possesses a frequency band structure exhibiting two characteristic features as can be observed in Fig. 3: two dispersive bands, which form a Dirac cone at the corners of the first Brillouin zone (for $\kappa_{x}=\kappa_{y}=\pm\pi$), and a flat band crossing the Dirac points. It is well-established that Dirac cones give rise to peculiar topological properties Weeks2010 and unusual behavior in general, such as effectively massless fermions, etc. Note that the flat-band frequency $\Omega_{FB}$ is equal to the resonance frequency of individual SQUIDs $\Omega_{SQ}$, i.e., $\Omega_{FB}=\Omega_{SQ}$. The maximum and minimum frequencies of the spectrum are obtained from Eq. (13) for $\kappa_{x}=\kappa_{y}=0$ and $\kappa_{x}=\kappa_{y}=\pi/2$, respectively, as

[TABLE]

Since $|\lambda_{x}|,|\lambda_{y}|\ll 1$, the bandwidth of the spectrum is approximatelly $\Delta\Omega\simeq 2\sqrt{\lambda_{x}^{2}+\lambda_{y}^{2}}/\Omega_{SQ}$. For example, for the parameters of Fig. 3 we have $\Omega_{min}\simeq 1.343$, $\Omega_{max}\simeq 1.384$, and $\Delta\Omega\simeq 0.04$. We also note that the flat band is an intrinsic property of this lattice in the nearest-neighbor coupling limit and thus it is not destroyed by any anisotropy (i.e., when $\lambda_{x}\neq\lambda_{y}$).

The dependence of the extremal frequencies $\Omega_{min,max}$ and the flat-band frequency $\Omega_{FB}$ on the parameters $\beta_{L}$ and $\lambda_{x},\lambda_{y}$ is shown in Fig. 4. In Fig. 4a all curves increase linearly with increasing $\beta_{L}$ while the bandwidth remains practically constant. In Fig. 4b the bandwidth increases with increasing $\lambda_{x}=\lambda_{y}$ while $\Omega_{FB}$ remains the same.

IV Numerics and Localization Measures

The dynamic equations for the fluxes through the SQUIDs, Eqs. (II), without losses and external forcing ($\gamma=0$ and $\phi_{ext}=0$), can be derived as the Hamilton’s equations from the Hamiltonian function

[TABLE]

where the energy (Hamiltonian) density $H_{n,m}$, defined as the energy per unit cell, is given by

[TABLE]

where $q_{n,m}^{k}=\frac{d\phi_{n,m}^{k}}{d\tau}$ is the normalized instantaneous voltage across the Josephson junction of the SQUID in the $(n,m)$th unit cell of kind $k$. Both $H$ and $H_{n,m}$ are normalized to the Josephson energy, $E_{J}$. The total energy $H$, given by Eqs. (15) and (IV), remains constant in time.

Eqs. (II) with $\gamma=0$ and $\Phi_{ext}=0$ are integrated in time using the second order symplectic Störmer-Verlet scheme Hairer2003 , which preserves the total energy $H$ to a prescribed accuracy which is a function of the time-step $h$. In the flux - voltage variables, that scheme reads Hairer2003

[TABLE]

where $\vec{\phi}^{k}$ and $\vec{q}^{k}$ are $N-$dimensional vectors ($N=N_{x}N_{y}$) containing the fluxes and the voltages for the SQUIDs of kind $k$ ($k=A,B,C$), and $H_{\vec{\phi}^{k}}\equiv\nabla_{\vec{\phi}^{k}}H$ denotes the column vector of partial derivatives of the Hamiltonian with respect to $\vec{\phi}^{k}$, i.e.,

[TABLE]

Periodic boundary conditions are used throughout, while the SLiMM is initialized at $\tau=0$ with a single-site excitation of amplitude $A_{m}$. The excited SQUID is either of kind $A$, $B$, or $C$. For isotropic coupling between SQUIDs, i.e., for $\lambda_{x}=\lambda_{y}$, a single-site excitation of either a $B$ or a $C$ SQUID provides identical results due to symmetry.

For the identification of the localized states that may be formed either due to the flat band or the nonlinearity, and the quantification of their degree of localization, two statistical measures will be used; the energetic participation ratio $P_{e}$ and the two-dimensional second moment $M_{2}$, which are given, respectively, by deMoura2003 ; Laptyeva2012 ; Mulansky2012

[TABLE]

and

[TABLE]

where $\epsilon_{n,m}=H_{n,m}/H$ is the normalized energy density, and $\bar{x}$, $\bar{y}$ are the coordinates of the ”center of energy”

[TABLE]

Note that $P_{e}$ measures roughly the number of excited cells in the system; its values range from $P_{e}=1$ (strong localization, all the energy in a single cell) to $P_{e}=N$, with $N=N_{x}N_{y}$ (equipartition of the energy over the $N$ SQUIDs). That measure has been also used to quantify the degree of diffraction in Kagomé photonic lattices Vicencio2014 . The second moment $M_{2}$ quantifies the squared width of the state, hence, its spreading.

Eqs. (II) with $\gamma=0$ and $\phi_{ext}=0$, implemented with periodic boundary conditions are initialized with single-site excitations of the form

[TABLE]

where $A_{m}$ is the amplitude of the initial excitation, and $k=A$, $B$ or $C$. The excited SQUID belongs to the unit cell with $n=n_{e},m=m_{e}$, with $n_{e}=N_{x}/2$ and $m_{e}=N_{y}/2$. The SLiMM is initialized with $A_{m}$ spanning several orders of magnitude, and for each $A_{m}$ several quantities such as the energy, the localization measures, and the ratio $r=|H(\tau)-H(0)|/H(0)$ are monitored during temporal evolution. Typically, a time-step $h=T_{SQ}/1000$, where $T_{SQ}=2\pi/\Omega_{SQ}$, is used in the simulations. However, it has been checked that smaller time-steps provide practically identical results. It has been also checked that in all runs the ratio $r$ remains less than $5\times 10^{-6}$ for the time step $h$ above.

V Flat-Band and Nonlinear Localization

The typical time-dependence of $P_{e}$ and $M_{2}$ when an edge SQUID (i.e., a SQUID $C$) is initialy excited with amplitude $A_{m}$ is shown in Figs. 5a and 5b, respectively, for $A_{m}=0.001$ (black), $0.01$ (red), $0.1$ (green), and $1$ (blue). Note that the curves for $A_{m}=0.001$ and $0.01$ almost coincide; for lower initial amplitudes the results are practically identical to those obtained for $A_{m}=0.001$. For such low initial amplitudes the SLiMM remains in the (almost) linear regime, in which localized flat-band states are expected to be observed. Indeed, as can be seen in Fig. 5a, as well as in the inset for $A_{m}=0.01$, $P_{e}$ has a running average over $5000~{}T_{SQ}$ time units which is about eleven ($P_{e}\simeq 11$, inset) indicating substantial localization. The existence of a localized state is advocated in this case by the corresponding second moment $M_{2}$, which running average over $5000~{}T_{SQ}$ time units (yellow curve) attains a constant value for relatively long integration times ($M_{2}\simeq 22$). The constancy of $M_{2}$ is interpreted as the termination of the energy spreading away from the site on which it was initially localized.

For $A_{m}=0.1$, the inspection of the corresponding (green) curve and its running average (maroon) reveals a dramatic change in the behavior of $P_{e}(\tau)$; the value of the latter increases more or less linearly with increasing $\tau$ until it saturates at a rather high value around $P_{e}\sim 140$. Note however the plateaus in the running average curve which indicate that the SLiMM passes through several metastable states until it reaches the steady one. The second moment $M_{2}$ in this case oscillates around $43$. Finally, for $A_{m}=1$ significant nonlinear effects come into play that favor strong localization with $P_{e}\sim 1$; thus, all the energy initially provided to the system at a single site, it practically remains there! This is actually the reason why the value of $M_{2}$ remains for all times close to zero ($M_{2}\simeq 0.1$, there is no spreading of energy whatsoever). Clearly, three different regimes can be identified; the (almost) linear regime, in which flat-band localization is possible, the intermediate regime, in which no localization is observed and the initial energy is eventually spread (in time-scales longer than those shown here) over the whole lattice, and the nonlinear regime in which localization in the form of intrinsically localized modes or discrete breathers, is observed. The size of fluctuations, e.g., in the curves for $P_{e}$, depends on that regime which in turn is determined by the initial condition (excitation); thus, fluctuations are weak in the linear, strong in the intermediate, and vanishing in the nonlinear regime.

The corresponding time-dependence of $P_{e}$ and $M_{2}$ when a corner SQUID (i.e., a SQUID $A$) is initialy excited with amplitude $A_{m}$ is shown in Figs. 6a and 6b, respectively. In this case, there is no localization in the linear and the intermediate regimes, i.e., for $A_{m}=0.001$ (black), $0.01$ (red), and $0.1$ (green), as can be inferred by the large values of $P_{e}$ whose running average over $5000~{}T_{SQ}$ time units is about $140$ ($P_{e}\sim 140$). At the initial stage of time integration which is not visible on the scale of the temporal axis of Fig. 6, both $P_{e}$ and $M_{2}$ have low values; however, within a few thousands time units they gradually grow to their high values. Note that the average of the curves for $M_{2}$ ($\sim 43$) is very close to that of the average of the corresponding curve for $A_{m}=0.1$ in Fig. 5b (intermediate regime). For high initial amplitude ($A_{m}=1$), however, strong localization due to nonlinearity is again observed. For such high values of $A_{m}$ the localized state which is generated either by initially exciting an $A$ or a $C$ SQUID does not reveal any significant difference. When there is no localization, the fluctuations of both $P_{e}$ and $M_{2}$ are again very strong. The results presented in Figs. 5 and 6 have been obtained for $\lambda_{x}=\lambda_{y}$, i.e., in the case of an isotropic Lieb lattice in the nearest-neighbor approximation. In this case, for single-site initial excitations of either a $C$ SQUID or a $B$ SQUID (i.e., of edge SQUIDs), the results are practically identical.

The above scenario is confirmed by inspecting the corresponding energy density plots, i.e., the plots of the energy density $E_{n,m}=H_{n,m}$ on the $n-m$ plane, which are shown in Fig. 7. In Fig. 7a, for $A_{m}=0.001$, the energy density $E_{n,m}$ is clearly localized, although not on only one unit cell; the energetic participation ratio is in this case $P_{e}\simeq 10.5$. A similar pattern is obtained for $A_{m}=0.01$, as shown in Fig. 7b, in which the maximum of the energy density is approximately two orders of magnitude larger than that in Fig. 7a. In Fig. 7c, there is clearly no localization, as it can also be inferred by the large participation ratio $P_{e}\simeq 140$. In Fig. 7d, in which the localization is due to the nonlinearity, the energy is almost completely localized, and $P_{e}\simeq 1$.

It should be noted that there are particular types of modes which cannot be efficiently excited in the SLiMM with initial single-site excitations used here. As an example, consider the application of the constraint $\phi_{n,m}^{A}=0$ for all the SQUIDs of kind $A$. That case has been also considered in a rhombic (quasi-) one-dimensional system with three waveguids per unit cell, whose coupling functions are the same with those of the equations for the SLiMM Maimistov2017 . By setting $\phi_{n,m}^{A}=0$ for all $n$ and $m$ and $\phi_{n,m}^{B}=\phi_{n,m}^{C}=\delta_{n,n_{1}}\delta_{m,m_{1}}$, with $n_{1}$ and $m_{1}$ integers, we get from the first of Eqs. (8) or the first of Eqs. (II) with $\gamma=0$ and $\phi_{ext}=0$, that $\lambda_{x}\phi_{n_{1},m_{1}}^{B}=-\lambda_{y}\phi_{n_{1},m_{1}}^{C}$ ($\phi_{n,m}^{B}=\phi_{n,m}^{C}=0$ for $n\neq n_{1}$ and $m\neq m_{1}$). That particular solution for the SLiMM system (either the linearized one or not) certainly cannot be obtained using single-site initial excitations.

In order to roughly determine the boundaries between the linear, intermediate, and nonlinear regimes, the averages of several quantities over the steady-state integration time $\tau_{int}$ were calculated for a wide range of initial excitation amplitudes $A_{m}=A_{m,i}$. An edge ($C$) SQUID is initially excited with amplitude $A_{m,i}$ and Eqs. (II) with $\gamma=0$ and $\phi_{ext}=0$ are integrated in time for $\tau_{int}=10^{5}~{}T_{SQ}$ time units, to allow for transients to die out (the obtained results are discarded) and the steady state to be reached. Then, in the steady state, the equations are integrated in time for $\tau_{int}$ more time units , and the energetic participation ratio averaged over $\tau_{int}$ is calculated. At the end of the steady-state integration time, the amplitude of the flux of the excited SQUID, $A_{m,c}$, and the oscillation frequency of the flux through the loop of the excited SQUID, $\Omega_{osc}$, are also calculated. The same calculations are performed for an initially excited corner SQUID $A$, and the results for both cases are shown in Fig. 8. In Fig. 8a, the calculated amplitude $A_{m,c}$ of the flux $\phi_{n_{e},m_{e}}^{k}$ through the loop of the SQUID with $k=A$ (blue) and $k=C$ (green) is shown along with an enlargement for low $A_{m,i}$ (inset). As it can be observed, $A_{m,c}$ attains low values for low initial amplitudes $A_{m,i}<0.15$, while for $A_{m,i}>0.15$ the calculated amplitude $A_{m,c}$ increases linearly with increasing $A_{m,i}$, according to the approximate relation $A_{m,c}\simeq A_{m,i}$. The behavior for $A_{m,i}>0.15$ is a result of the strong localization taking place due to nonlinearities and it does not depend on which kind of SQUID (edge or corner) is initially excited. However, a closer look to the two curves for $A_{m,i}<0.15$, reveals significant differences, especially for $A_{m,i}<0.05$, which can be seen more clearly in the inset. In this regime the calculated amplitude $A_{m,c}$ for $k=C$ follows the relation $A_{m,c}\simeq A_{m,i}/2$, indicating localization due to the flat band. This conclusion is also supported by Figs. 8b and 8c. In Fig. 8b, the energetic participation ratio averaged over $\tau_{int}$, $<P_{e}>$, for low values of $A_{m,i}$ attains very different values depending on which kind of SQUID is initially excited ($A$ or $C$); specifically, while $<P_{e}>\sim 10.5$ for the SLiMM when a $C$ SQUID is initially excited (black), it is $<P_{e}>\sim 140$ when an $A$ SQUID is initially excited (red). That large difference between the values of $<P_{e}>$ is due to flat-band localization in the former case and delocalization in the latter case since no flat-band modes are excited. In the inset, it can be observed that $<P_{e}>$ for an initially excited $C$ SQUID starts increasing for $A_{m,i}>0.05$ indicating gradual degradation of flat-band localization and meets the $<P_{e}>$ curve for an initially excited $A$ SQUID at $A_{m,i}\sim 0.1$. In Fig. 8c, for $A_{m,i}<0.15$, the oscillation frequency of the flux through the initially excited SQUID $\Omega_{osc}$ (either $A$ or $C$), has a value around that of the linear resonance frequency of a single SQUID, $\Omega_{SQ}$ ($\Omega_{SQ}\simeq 1.364$ for the parameters of Fig. 8). As it can be seen in the inset, when a $C$ SQUID is initially excited (violet), then up to high accuracy $\Omega_{osc}=\Omega_{SQ}$ for initial amplitudes up to $A_{m,i}\sim 0.075$. However, when an $A$ SQUID is initially excited (turquise), the frequency $\Omega_{osc}$ jumps slightly above and below $\Omega_{SQ}$ irregularly, but it remains within the bandwidth of the linear frequency spectrum. For $A_{m,i}>0.15$, the frequency $\Omega_{osc}$ decreases with increasing $A_{m,i}$, although it starts increasing again with increasing $A_{m,i}$ at $A_{m,i}\sim 0.8$. In this regime, nonlinear localized modes of the breather type are formed, which frequency lies outside the linear frequency spectrum and depends on its amplitude, as it should be. From this figure it can thus be inferred that flat-band localization occurs for initial amplitudes up to $A_{m,i}\simeq 0.05$ (linear regime), while delocalization occurs in the interval $0.05<A_{m,i}<0.15$ (intermediate regime). For larger $A_{m,i}$, strong nonlinear localization occurs (nonlinear regime). This rough estimation for the boundaries between the three regimes is of course parameter dependent. Remarkably, flat-band localization occurs only when an edge SQUID ($B$ or $C$) is initially excited. The excitation of a corner ($A$) SQUID does not lead to excitation of flat-band modes and thus such a localized initial state rapidly delocalizes. On the other hand, the observed flat-band localization is not very strong as compared to the nonlinear localization. This is probably due to the fact that a single-site excitation of a $B$ or $C$ SQUID does not correspond to an exact localized flat-band eigenmode.

In the case of anisotropic coupling, i.e., for $\lambda_{x}\neq\lambda_{y}$, single-site excitations of $B$ and $C$ SQUIDs give different results as expected due to the lowering of symmetry Guzman-Silva2014 . Typical curves for the amplitude of the flux $\phi_{n_{e},m_{e}}^{k}$ of the excited SQUID $A_{m,c}$ ($k=A$, $B$, and $C$), the energetic participation ratio averaged over the steady-state integration time $<P_{e}>$, and the oscillation frequency $\Omega_{osc}$ of the flux through the excited SQUID, are shown in Fig. 9 for anisotropic nearest-neighbor coupling, $\lambda_{y}=1.5\,\lambda_{x}=-0.03$, as a function of relatively low initial excitation amplitudes $A_{m,i}$ for which flat-band localization is expected. As can be observed in Fig. 9a, flat-band localization occurs when either of the edge SQUIDs are excited with $A_{m,i}<0.04$. For larger values of $A_{m,i}$, localization starts degrading as it can be confirmed from the corresponding curves of $<P_{e}>$ in Fig. 9b. Here, it is also apparent that initial excitations of $B$ and $C$ SQUIDs do not lead to a state with the same degree of localization; indeed, $<P_{e}>$ is respectively $\sim 5$ and $\sim 25$ (while the corresponding $<P_{e}>$ for initial excitations of an $A$ SQUID is about $160)$. The corresponding oscillation frequencies of the fluxes in the case of $B$ or $C$ SQUID initial excitations are practically equal to that of the single SQUID resonance $\Omega_{SQ}\simeq 1.364$. For initial excitations of an $A$ SQUID, the oscillation frequency is very close to that of either the upper or the lower boundary of the linear frequency spectrum.

VI Conclusions

The dynamic equations for the fluxes threading the SQUID loops of a driven-dissipative SLiMM have been derived, along with the corresponding linear frequency spectrum. The Lieb lattice geometry results in a spectrum with two dispersive bands, which form a Dirac cone at the corners of the first Brillouin zone, and a flat band crossing those Dirac points. The localization properties of Hamiltonian SLiMMs, i.e., those without dissipation and driving terms, have been determined through numerical simulations for single-site initial excitations of varying amplitude. Flat-band localization, i.e., the emergence of localized flat-band states, is observed when an edge ($B$ or $C$) SQUID of the unit cell of the SLiMM is initially excited with low amplitude. To the contrary, no such states are generated when a corner ($A$) SQUID of the unit cell of the SLiMM is initially excited with low amplitude. These results are compatible with the experiments on photonic Lieb lattices Mukherjee2015a . For sufficiently high amplitude of the initial excitation of either a corner or an edge SQUID, localization due to nonlinearities in the form of discrete breathers is observed. The linear (low amplitude initial excitations) and the nonlinear regimes (high amplitude initial excitations), in which flat-band localized states and discrete breathers, respectively, can be generated, are separated by an intermediate regime in which neither type of localization is observed. This dynamic behavior is quite different from that observed in, e.g., two-dimensional Kagomé lattices in which families of nonlinear localized modes in the form of discrete solitons or discrete breathers may bifurcate from localized linear modes of the flat band Vicencio2013 ; Johansson2015 . Here, relatively high-amplitude initial excitations ($A_{m,i}>0.05$) excite nonlinear effects in the SQUIDs which destroy the flatness of the flat-band which has been obtained in the linear limit. At the same time, however, these nonlinear effects are not strong enough to help the initial excitation to remain localized (self-trapped); that occurs only when the amplitude of the initial excitation exceeds a particular, parameter-dependent threshold ($A_{m,i}\simeq 0.15$ for the parameters of Fig. 8).

ACKNOWLEDGMENT

This work is partially supported by the Ministry of Education and Science of the Russian Federation in the framework of the Increase Competitiveness Program of NUST ”MISiS” (No. K2-2015-007) and by the Ministry of Education and Science of the Republic of Kazakhstan (Contract $\#$ 339/76 –2015). NL gratefully acknowledges the Laboratory for Superconducting Metamaterials, NUST ”MISiS” for its warm hospitality during visits.