Floquet bound states around defects and adatoms in graphene

Recent studies have focused on laser-induced gaps in graphene which have been shown to have a topological origin, thereby hosting robust states at the sample edges. While the focus has remained mainly on these topological chiral edge states, the Floquet bound states around defects lack a detailed study. In this paper we present such a study covering large defects of different shape and also vacancy-like defects and adatoms at the dynamical gap at $\hbar\Omega/2$ ($\hbar\Omega$ being the photon energy). Our results, based on analytical calculations as well as numerics for full tight-binding models, show that the bound states are chiral and appear in a number which grows with the defect size. Furthermore, while the bound states exist regardless the type of the defect's edge termination (zigzag, armchair, mixed), the spectrum is strongly dependent on it. In the case of top adatoms, the bound states quasi-energies depend on the adatoms energy. The appearance of such bound states might open the door to the presence of topological effects on the bulk transport properties of dirty graphene.


I. INTRODUCTION
Driving a material out of equilibrium offers interesting paths to alter and tune its electrical response. A prominent example is the generation of light-induced topological properties, 1-3 e.g. illuminating a material like graphene to transform it in a Floquet topological insulator (FTI). Very much as ordinary topological insulators (TI), [4][5][6][7] FTIs have a gap in their bulk (quasi-) energy spectrum-being then a bulk insulatorand their Floquet-Bloch bands are characterized by non-trivial topological invariants. 3,8,9 In addition, and despite some important differences with TIs, 8,10 FTIs show a bulk-boundary correspondence and hence host chiral/helical states at the sample boundaries.
The emergence of such non-equilibrium properties has been intensively investigated in recent years in a variety of systems including graphene [11][12][13][14][15][16][17][18][19] and other 2D materials, 20,21 normal insulators, 2,22 coupled Rashba wires, 23 photonic crystals, 24 cold atoms in optical lattices, [25][26][27][28][29][30][31] topological insulators, [32][33][34][35][36] and also classical systems. 37 The research interest has focused in many different aspects of the problem such as the characterization of the edge states, 16,17 different signatures in magnetization and tunneling, 38,39 the proper invariants entering the bulk-boundary correspondence, 8,10,19,40 their statistical properties, 41,42 the role of interactions and dissipation 41,[43][44][45][46] and the associated two-terminal 47,48 and multiterminal (Hall) conductance both in the scattering 49 and decoherent regimes. 45 So far, however, the experimental confirmation of the presence of such edge states has only been achieved in photonic crystals. 24 Nonetheless, in condensed matter systems the Floquet induced gaps have already been observed at the surface of a topological insulator (Bi 2 Se 3 ) by using time and angle resolved photoemission spectroscopy (tr-ARPES). 33 More recently, effective Floquet Hamiltonians were realized in cold matter systems. 50 Despite the intense research on FTIs, most of the studies address pristine samples. Besides occurring naturally in any sample, defects will also host Floquet bound states when the sample is illuminated. If the defects are extended, the presence of the associated Floquet bound states might allow for new experiments probing them. This motivates our present study. Specifically, taking laser-illuminated graphene as a paradigmatic example of a FTI, we study Floquet bound states around defects in the bulk of a sample. We show that chiral states circulate around holes or multi-vacancy defects of different shapes and lattice terminations (zigzag, armchair or mixed) like the ones showed in Fig. 1. The properties of these states (quasi-energies and their scaling with the system parameters, associated probability currents, etc.) are characterized using both numerical simulations, by means of a tightbinding model, and analytical approaches, by solving the appropriate low energy Dirac Hamiltonian in a reduced Floquet space. Quite interestingly, these bound states persist even in the limit of a single vacancy defect. Furthermore, bound states are found around adatoms that sit on top of a C atom (like H or F, for instance).
While the presence of Floquet bound states around vacancy-like defects or adatoms might jeopardize the exper- imental observation of laser-induced gaps, they could, on the other hand, also open the route towards the observation of interesting topological transport phenomena in dirty bulk samples by changing localization or percolation properties, for instance.
The rest of the paper is organized as follows. First, we introduce our low energy model and the associated analytical Floquet solutions (Sec. II). Several particular cases are presented in section III, namely, large holes with zigzag or armchair edge terminations, as well as defects consisting of regions with a staggered potential. The chiral nature of the currents associated to the bound states is discussed in Sec. IV. In Sec. V we compare our solutions with numerical calculations on a tight-binding model. The case of point like defects such as vacancies or adatoms is presented in Sec. VI. We finally conclude in Sec. VII.

II. THE LOW ENERGY MODEL AND THE FLOQUET SOLUTION
Let us consider an irradiated graphene sample with a single defect. Since the bound states we want to describe are topological in origin, 16,17,19 the specific form or nature of the defect (see Fig. 1) is irrelevant for probing their existencethough the details of the quasi-energy spectrum and the particular form of the wave-functions will depend on it. To simplify the discussion we will start by assuming that the defect potential does not mix the different graphene valleys (Dirac cones)-this assumption will be relaxed when discussing particular examples. Hence, the low energy behavior around both cones can be described by a Hamiltonian given bŷ if we use the isotropic representation where the K and K ′ cones are described by the wave-functions Here v F ≃ 10 6 m/s denotes the Fermi velocity, σ = (σ x , σ y ) represents the Pauli matrices describing the pseudo-spin degree of freedom (sites A and B of the honeycomb lattice), e is the absolute value of the electron charge, c is the speed of light and A(t) = Re A 0 e iΩt the vector potential of the electromagnetic field (a plane wave incident perpendicularly to the graphene sheet). The associated electric field is then It is important to emphasize that while we will refer to graphene from hereon, our results apply to any massless Dirac fermion system described by Eq. (1).
Since for solving the time-dependent Schrödinger equation we will take advantage of the Floquet formalism 51,52 used to deal with time dependent periodic Hamiltonians, it is instructive to briefly introduce its basic ideas (for a more extensive general reviews we refer to Refs. [53] and [54]). Floquet theorem guarantees the existence of a set of solutions of the form ψ α (t)⟩ = exp(−iε α t ̵ h) φ α (t)⟩ where φ α (t)⟩ has the same time-periodicity as the Hamiltonian, φ α (t + T )⟩ = φ α (t)⟩ with T = 2π Ω. 51,53 The Floquet states φ α ⟩ are the solutions of the equationĤ whereĤ F =Ĥ − i ̵ h∂ t is the Floquet Hamiltonian and ε α the quasi-energy. Using the fact that the Floquet eigenfunctions are periodic in time, it is customary to introduce an extended R ⊗ T space (the Floquet or Sambe space 52 ), where R is the usual Hilbert space and T is the space of periodic functions with period T . A convenient basis of R ⊗ T can be built from the product of an arbitrary basis of R (the eigenfunctions a n ⟩ of the time-independent part of the Hamiltonian, for instance) and the set of orthonormal functions e imΩt , with m = 0, ±1, ±2, ... that span T . Then, or, in a vector notation in R ⊗ T , Here, u α m ⟩ = ∑ n B α mn a n ⟩ are linear combinations of the basis states of R. Written in this basis,Ĥ F is a time-independent infinite matrix operator with Floquet replicas shifted by a diagonal term m ̵ hΩ and coupled by the radiation field with the condition, for pure harmonic potentials, that ∆m = ±1.
In the absence of any defect, the Floquet spectrum presents dynamical gaps at different quasi-energies. 1,17,19 Here, we will focus on the gap, of order η ̵ hΩ, that appears at ε ∼ ̵ hΩ 2 and look for bound states inside it. Since we will only consider the limit η = v F eA 0 c ̵ hΩ ≪ 1, it is sufficient to restrict the Floquet Hamiltonian to the m = 0 and m = 1 subspaces (or replicas) for the analytical calculations-the numerical results can retain a larger number (N FR ) of replicas if necessary. As discussed in Refs. [17] and [19], this restriction is enough to get the main features of the energy dispersion and the Floquet states when η ≪ 1.
The reduced Floquet Hamiltonian describing states near ε ∼ ̵ hΩ 2 then corresponds tõ It is straightforward to see thatH F φ(r) = εφ(r) implies that and hence only two functions, u 0A (r) and u 1B (r), have to be found. These functions satisfy where p 2 = p + p − = p − p + . Because we are interested in describing the effect of a defect-which breaks the translational invariance of the systems-, it is useful to change at this point to a polar coordinate system, r and ϕ, centered at it. In terms of these variables we have, Similarly, as in the case of local defects in ordinary TI, 55,56 the solutions of Eq. (8) can be written as u 1B (r) = e ilϕ f (k 0 r) and u 0A (r) = e ilϕ g(k 0 r) with l an integer number. This follows from the fact that [H F , L] = 0, where and L φ(r) = ̵ hl φ(r), where φ(r) is given by Eq. (6). In order to proceed further we define the adimensional parameters With this notation, the equations for f (ξ) and g(ξ) become For quasi-energies inside the bulk dynamical gap, the wavefunction must decay far from the defect. Hence, let us look for a solution of the form f (ξ) = c K l (λξ) and g(ξ) = d K l (λξ), where K l (x) is the modified Bessel function of the 2nd kind that satisfy Introducing this into Eqs. (12) we arrive to the following condition for λ, and the relation The equation for λ has four solutions which are complex conjugate in pairs. The two physical solutions correspond to Re(λ) > 0 as this guarantees an exponential decay for large r.
Let us denote these two solutions as λ + and λ − = λ * + , The region where Re(λ) > 0 corresponds to µ < η 1 + η 2 , that is, inside the bulk dynamical gap, 17 ∆ = ̵ hΩη 1 + η 2 . The other components of the Floquet wavefunction can be readily obtained as which are straightforward to evaluate since ∂ ξ ∓ l ξ K l (λξ) = −λK l±1 (λξ). It is worth to point out that ⟨u 1 u 0 ⟩ = 0 so that φ(r, t) can be normalized for any time t in this approximation, 17 which allows to calculate not only time-averaged quantities but also their time dependence explicitly.
To proceed any further we need to specify the defect type, which allows the setting of the appropriate boundary conditions. In the following we present a detailed discussion for some particular but relevant cases.

III. BOUNDARY CONDITIONS
The boundary conditions (BC) must guarantee that the probability current perpendicular to the defect boundary cancels out. Here, we shall consider only three types of BCs that represent three generic cases and serve to illustrate the overall picture: the zigzag-like BC (ZZBC), the armchair-like BC (ABC) and the infinite mass BC (IMBC). 57 Since the BC needs to be satisfied at any time, in Floquet space the boundary condition must be imposed on each replica separately. Therefore, the boundary problem is analogous to the static one and we shall follow Refs. [58] and [59] and use a matrix M to introduce the appropriate relations between the components of the A and B sublattices and the two Dirac cones at the boundary for the three types of BCs. 59,60 An arbitrary BC can be written in the form where R(ϕ) defines the shape of the defect and the matrix M (in the isotropic representation) is given by Here σ refers to the sublattice pseudospin and τ to the valley (Dirac cones) isospin. The matrix M has all the information about the shape of the boundary via the unit vectorn. On the other hand, the nature of the honeycomb lattice's termination is related to the unit vectorν, that rules whether the two Dirac cones mix or not. Namely, for a defect with a straight boundary, 59 ZZBC →ν =ẑ ,n = ±ẑ ABC →ν ⋅ẑ = 0 ,n =ẑ ×n B (20) IMBC →ν =ẑ ,n =ẑ ×n B , wheren B is an unitary vector perpendicular to the defect boundary and pointing inwards. From the above expressions it is clear that while armchair BC mixes cones, zigzag and infinite mass BCs do not. In the following we shall be interested in the comparison between analytical and numerical results for simple geometries, and so we will restrict ourselves to handle only defects with regular polygonal shapes with N sides. The general form of M N for such cases is given in the Appendix A.
While for the honeycomb lattice, defects with well defined terminations can only have N = 3 or N = 6, it is useful to discuss the limiting case of a circular defect and then compare with the numerics. For the ABC and IMBC this corresponds to the limit N → ∞ while for the ZZBC care is needed to account for the change of the sublattice character of the edge atoms [n = ±ẑ depending on the sublattice].
A. Circular defect with "zigzag" boundary condition The ZZBC does not mix valleys. This is valid for arbitrary N , i.e. M N is diagonal in the isospin subspace. Moreover, it is also diagonal in the pseudospin subspace. However, it is possible, as in the hexagonal geometry, that different sides of the polygon terminate in sites corresponding to different sublattices. This is represented by then = ±ẑ in Eq. (20), where the sign changes from side to side, thereby making it cumbersome to handle analytically. Hence, for the sake of simplicity, we will consider a 'fictitious' case where the ± sign is ignored and later compare with the exact numerical calculation. Hereon we will refer to it as the circular-ZZBC (cZZBC).This will help us to better grasp some aspects of the problem.
For a circular defect (of radius R) the BC implies, say, that u 1B ( r = R) = 0 and u 0B ( r = R) = 0-this corresponds to a honeycomb lattice that ends on A sites. To satisfy it we need to combine the two independent bulk solutions discussed in Section II. That is, where we have kept the previous notation. Then we have that with ξ 0 = k 0 R. This leads to the following relations between coefficients: c + = c − and d + = d − . By introducing them back into Eqs. (12) we obtain, for the K cone, the following equation for the quasi-energy (µ) with The solutions (µ l ) to this equation form a discrete set of quasienergies inside the bulk dynamical gap. Figure 2 shows them as a function of ξ 0 (throughout this work, we shall use η = . Notice that the symmetry between l > 0 and l < 0 is broken by the radiation field. The symmetry of the Floquet spectrum around the center of the gap (µ = 0) is recovered when the complementary valley (K ′ cone) is considered. For that, we recall that the solutions for the K ′ cone can be obtained by relabeling the Floquet wavefunction as . This results in an additional set of quasi-energies that can be obtained from the condition It can be shown that the latter set of quasi-energies can be obtained from Eq. (23) by exchanging (l, µ) → (−l, −µ), which is precisely what is needed to recover the symmetry around µ = 0. It is interesting to consider, for a fixed l, the limit of very large radii, ξ 0 ≫ ξ d = k 0 ̵ hv F ∆ = 1 + η 2 2η and approximate K l (λξ 0 ) by its asymptotic expansion. By doing so, Eqs. (23) and Eq. (25) leads to respectively. This result can be understood in terms of the quasi-energy dispersion of the edge states in irradiated semiinfinite graphene sheets with a zigzag termination. 17 In that case, it was shown that, close to the center of the gap, the quasi-energy dispersion can be approximated by ε k = ̵ hΩ 2 ± ̵ hΩη 2 2 + ̵ hv F ηk. Our result for µ l is then reflecting the fact that the wavevector k along the defect's edge must be quantized, It is worth mentioning that in this large radii limit the Floquet states have roughly the same weigth on the two Floquet replicas.

B. Infinite mass boundary condition
The IMBC was introduced by Berry and Mondragon in Ref. [57] to study confined Dirac particles ('neutrino billiards'). It corresponds to add a mass term to the Dirac equation only in a given region of space (in our case the defect) and take the limit of such a mass going to infinity. While this could be thought as a local staggered potential in the honeycomb lattice, it must be kept in mind that this is only the case for a staggered po-tential much smaller than the bandwidth-this is so because if the staggered potential is too large it behaves like an effective hole (introducing inter-valley scattering depending on the geometry of the defect). The latter limit was not a problem in Ref. [57] , because they only considered a single unbound massless Dirac particle.
Since the IMBC does not mix valleys either, we can treat again both Dirac cones separately. We start by using the circular geometry, which corresponds to the N → ∞ limit of M N . For the IMBC M ∞ is not longer diagonal in the pseudospin subspace and thus the A and B components of the wavefunction are not independent any more. In fact, Eq. (18) requires that 57 for the K and K ′ cone, respectively, where j = 0, 1 is the Following the same procedure as in the previous section, and using the same notation, these conditions imply that while the equation for the quasi-energies is given by Here the (-) and (+) signs correspond to K and K ′ cone, respectively. It can be shown that the above expression remains invariant under the change (µ, l) → (−µ, −l) for each cone separately and, therefore, unlike the cZZBC, the Floquet spectrum for the IMBC is symmetric around µ = 0 for each cone. Using this symmetry of Eq. (30) it is straightforward to verify that there is no solution for l = 0 (that necessarily corresponds to µ = 0). The IMBC Floquet spectrum is shown in Fig. 3 as a function of ξ 0 . Note that the two cones have a completely different spectrum. This could be anticipated from the fact that the presence of both the staggered potential and the radiation field breaks the valley symmetry (cf. Fig. 5 below)-it is worth mentioning that the bulk Floquet gap at k = 0 can even present a topological phase transition depending on the relative magnitude of the mass term and the radiation field. 61 When defects are made of regular polygons, i.e. with finite N , the M N matrix acquire a non-trivial structure as a function of ϕ. Thus, the states whose quantum numbers l differ in N are coupled, thereby leading to avoided crossings. The equations for this case are rather cumbersome (some of them are presented in the appendix) but can be solved in a perturbative fashion. Some examples are presented in Sec V in comparison with the numerical solutions of the tight-binding model.

C. Armchair boundary condition
The ACB is analog to the IMBC in the pseudospin subspace, leading to similar quasi-energy spectra. The difference between both boundary conditions rely on the isospin subspace: while ACB mixes cones, IMBC does not. Thus, ACB exhibits additional avoided crossings between modes belonging to different cones (see numerical results in Sec. V). Because cones are mixed, they both need to be treated together and hence the dimension of the Floquet space is doubled. The analytical procedure is similar to the one presented for the other BCs, whose details are beyond the scope of the present work. We will then limit, for this case, to discuss the numerical results in in Sec. V.

IV. PROBABILITY CURRENT DENSITY: CHIRAL CURRENT
So far we have mainly analyzed the spectrum of the Floquet bound states inside the dynamical gap (around ̵ hΩ 2) for a circular defect. Now we focus on their chiral nature. The velocity operator is given byv = v F σ and hence the time averaged (over one period) probability current density is where ⟨σ α ⟩ j = {u * jA,l (r), u * jB,l (r)}σ α {u jA,l (r), u jB,l (r)} T , j = 0, 1 is the same as earlier, σ r = σ ⋅r and σ ϕ = σ ⋅φ. Using the solutions founded in the previous section, it can be readily shown that Since λ + = λ * − , one can easily check that Im f l (ξ)f * l (ξ) = Im (g l (ξ)g * l (ξ)) = 0 so that the radial component of the current density vanishes, as expected. Therefore, we have Figure 4 shows the spatial dependence of both the probability and the current density for the K and K ′ cones and for the two different boundary conditions analyzed in Sec. III. The curves correspond to a defect of R = 30 a cc , i.e., ξ 0 = 1 with the parameters used throughout this work. We have only retained the Floquet wavefunctions with l = 0, 1, 2, whose corresponding quasi-energies can be seen from Fig. 2 and Fig. 3 for ξ 0 = 1. Due to the oscillating nature of the Floquet wavefunctions both probability density functions and current densities show relative maxima and minima (with the same or different signs in the case of current densities) as a function of ξ. Nevertheless, all of them decay exponentially away from the edge of the defect. This is more evident for the Floquet wavefunctions whose quasi-energies are close to the middle of the dynamical gap as in that case the decay length is shorter. For quasi-energies close to the edges of the dynamical gap, the decay length becomes larger and larger and the ξ −1 2 power law decay, characteristic of the K l Bessel functions with purely imaginary argument becomes apparent. In these latter cases, however, the current amplitude becomes several orders of magnitude smaller than in the formers (see Fig. 6). For the cZZBC, Fig. 4 shows the equivalent role that play the K and K ′ cones under the change l ↔ −l, as it was explained before in Sec. III A. Unlike the cZZBC, for the IMBC the K and K ′ cones are inequivalent. In this case, as discussed in Sec. III B, the change l ↔ −l lead to the same probability and current densities for each cone separately.
The lack of equivalence between the K and K ′ cones for defects with IMBC is also present in systems other than circular defects. For illustrating purposes, Fig. 5 shows the kdependent local density of states (LDOS) for a nanoribbon with both cZZBC and IMBC, projected on the m = 0 Floquet replica. Notice that, unlike the cZZBC, the IMBC presents an asymmetry (at each edge) with respect to the middle of the dynamical gap. The symmetry is broken by the presence of the mass term at the edges and it is only globally recovered when both edges are considered-this is so because for zigzag nanoribbons, as considered here, the atoms at the two edges belong to different sublattices.
Even when the current density oscillates as it decays away from the defect, the total current (current densities integrated  on r) for cZZBC has the same sign for all the bound states. This is the signature of the chirality of the Floquet states and their signs only depends on the sign of the helicity of the circularly polarized radiation field. Figure 6 shows the total currents for both cZZBC and IMBC as a function of the quantum number l for defects with ξ 0 = 1, 5, 10, 20. Unlike the cZZBC, the IMBC only presents chiral Floquet states for the K cone. Analogously, Fig. 5 shows a similar behavior for the nanoribbon with IMBC: while the K cone presents two chiral states at each edge, K ′ cone has none.
Finally, it is interesting to analyze the value of the total current of a given bound state in the limit of a large defect. As discussed in Sec. III A for large R the quasi-energy dispersion can be related to the one corresponding to a nanoribbon as the boundary of the defect appears (locally) as a straight line (i.e. when the radius is much larger than the decay length). In that case the expected velocity for each bound states is . 17 The inset of the Fig. 6 shows the current in units of v F for Floquet states with l = 0 (red points) as a function of the size of the defects. The black dotted line represent the expected η (1 + η 2 )-this is also indicated in the main figure. Clearly, there is a good agreement with the expected value. A similar behavior is observed for states with different quantum number l as the size of the defect increases.

V. COMPARISON WITH THE TIGHT-BINDING MODEL
In this section, we calculate the quasi-energy spectra within the dynamical gap numerically as a function of the size and shape of the defect for all three types of boundary conditions mentioned before, ZZBC, ABC and IMBC, and compare with the analytical results when possible.
In order to describe the electronic structure of irradiated graphene sheets near the Fermi energy, we resort to the widely used tight-binding Hamiltonian, [62][63][64] which is written only in terms of p z orbitals with energies i for a given carbon atom located at site i and hopping matrix elements γ ij between nearest-neighbors carbon atoms. In second quantization notation, it results where the operator c † i (c i ) creates (annihilates) a p z -electron on site i. The effect of the laser is introduced through the time-dependent phase of the hopping matrix elements, 1,65,66 where Φ 0 is the magnetic flux quantum and γ 0 ∼ 2.7 eV. 67 By using Floquet theory 54,68,69 as described before one can compute the Floquet spectrum. Once again, one ends up with a time-independent problem in an expanded space. In this case one can picture it as tight-binding problem in a multichannel system where each channel represents the graphene sheet with different number of photons. 51,66,70 It is worth mentioning that in the tight-binding method the time dependent perturbation is never purely harmonic given the exponential dependence of Eq. (35) on the radiation field amplitude. Hence, there is a coupling among all the replicas 66 and not just those with ∆m = ±1. Nevertheless, for η ≪ 1, only the latter are relevant.
Because the problem in the Floquet space becomes time independent, one can use standard techniques to calculate the quasi-energy spectrum. In this case we used the Chebyshev's polynomials method 71 which provides an order N method of proven efficiency. 72 This allows us to tackle very large systems sizes so that our defect is far from the boundaries and can be considered as a 'bulk defect'. For simplicity we only retained two Floquet replicas just like its theoretical counterpart studied in Section II. This is a good approximation whenever η ≪ 1. The addition of more replicas would lead to the development of a hierarchy of bound states in a similar way as for edge states at the border of an irradiated graphene sample. 19 Defects were introduced in graphene by defining geometrical shapes-triangles, hexagons, and circles-and removing all atoms inside it (for the ZZBC and ABC) as well as any remaining dangling bonds. In the case of the IMBC, a staggered potential was introduced only inside the defect-i.e. we added on-site energies (±δ) whose signs depend on the sublattice index. In all calculations we used δ = γ 0 2, which is larger than ̵ hΩ 2 (taken to be ∼ γ 0 20) but not too large as to become equivalent to a hole (δ → ∞ is equivalent to a hole defect). Triangles and hexagons in arbitrary orientations lead to edges with mixed zigzag and armchair terminations. However, for specific orientations with respect to the C-C bonds, it is possible to construct defects with only one termination type-we will refer to them as zigzag/armchair triangular and hexagonal defects. Circles, of course, are always a mixture of different edge terminations and, as we will show, present some special features. In all cases, the numerical calculations were performed using graphene samples of 1000 × 1000 unit cells. Figures 7 and 8 show a color map of the Floquet local density of states (FLDOS) inside the bulk gap (projected onto a few sites around the defect boundary, and on the m = 0 replica) as a function of the size of the defect for hole and staggered potential defects, respectively. The shape of the defect is indicated in the figures. Left panels correspond to zigzag terminations and the right panels to the armchair ones. Dashed (black) lines correspond to the solutions obtained from the continuum model (see discussion below). It is apparent from the figures that discrete Floquet bound states do appear inside the dynamical gap. Interestingly, in most cases, the quasi-energy spectrum resemble the ones obtained with the analytical model proposed in Sec. II. This remains valid for the triangular shaped zigzag hole even when the analytical solution relies on the circular symmetry of the defects. It is worth mentioning that for a quantitative comparison an effective radius is needed. In these cases we used R = 1 (2π) ∫ radius R(ϕ), as discussed in Sec. III and the appendix. This avoided crossings occurs whenever the quantum numbers of the crossing levels, l and l ′ , differ in a multiple of the number of sides N . A few particular examples are indicated in the Fig. 8. (ii) the latter picture is very particular in the case of the zigzag triangular hole defect (top-left in the Fig. 7). On the one hand, the matrix M is independent of ϕ-note thatn =ẑ for any ϕ as the edge site always belong to the same sublattice and the direction ofν is fixed for each cone-and hence the only dependence on ϕ appears through the boundary radius R(ϕ). On the other hand, for each cone, the 'unperturbed' energy levels of the 'zigzag circle' are never degenerated, making the effect even weaker. As a result, the energy level are well described by assuming that there is no mixing between states with different quantum number l. Notice also there is no mixing between different cones or valleys.
(iii) the zigzag triangular defect with the staggered potential shows a shift in energy with respect to the IMBC solution. This is related to the sublattice imbalance of the edge sites and the fact that both sublattices have different energy inside the defect (staggered potential). This effect is not observed for the other geometries as they have balanced edges.
(iv) the armchair hexagonal hole defect shows two distinct contributions to the quasi-energy spectrum. The one shown in Fig. 7, that is very close to the analytical solution for IMBC [except for the anticrossings between energy levels belonging to different cones that are only present in the armchair case (black arrows)], and the one presented in Fig. 14  completely different pattern. The two cases differ in the way the atom chains that constitute each side match at the vertices.
(v) the zigzag hexagonal hole defect presents a rather complex spectrum quite different from the rest. This is related to the strong mixing between states with different l imposed by the BC that requires that alternating components of the wavefunction cancel in alternating sides. A precise description of this case is beyond the scope of the present work.
Finally, we show numerical results for circular defects in the Fig. 9. The top panel corresponds to a hole defect and the bottom one to the staggered potential defect. Clearly, the latter is very well described by the analytical solutions (dashed black lines). Notice that no avoided crossings (if they exist) are resolved in our numeric simulations, presumably because they are very small since the actual geometry of the defect is very close to a circle. The spectrum of the circular hole defect is, as in the zigzag hexagonal one, very complex. Here, however, a more regular pattern emerges for large R as the quasi-energy of the bound states are pretty much confined to regions delimited by the analytical solution of the zigzag circular defect (dashed lines).
One of the questions that remains is to what extent do these bound states survive in the limit of a vacancy defect or, more generally, in the case of adatoms. This is particularly important as the presence of bound states around such impurities might hinder the ability to resolve the laser-induced gaps in actual experiments or lead to percolating states in dirty samples.

VI. THE ADATOM AND VACANCY DEFECTS
The continuum model presented in Sec. II is not adequate for analyzing the vacancy limit. In fact, in the R → 0 limit for zigzag hole (the appropriate one for a vacancy defect) one finds that there are no solutions inside the gap. Of course, this is not the correct approach as one should introduce a spatial cutoff to account for the finite size of the defect. In this sense, a tight-binding model approach is more convenient and allows for its generalization to include the adatom case.
Since we focus on the bound states within the dynamical gap at ̵ hΩ 2, it is enough to consider, as before, only two Floquet replicas, m = 0 and m = 1. While for the numerical calculations we will use the real space version of the tightbinding Hamiltonian presented in the previous section, for the discussion of the main aspects of the problem it is better to use a k-space representation. Then, the Floquet Hamiltonian is written as Here a † mk and b † mk create an electron on the Floquet replica m on the Bloch state with momentum k on the sublattice A and B, respectively, φ k = ∑ δj e jk⋅δj , where {δ i } are the relative coordinates of the three nearest neighbors A sites of a given B site, t = γ 0 J 0 (z) , and A k = γ 0 J 1 (z) ∑ δj e ik⋅δj (δ jx − iδ jy ) a cc with J n (x) the n-th Bessel function of the first kind and z = 2πA 0 a cc Φ 0 . 66 We describe the adatom impurity with a single orbital of energy bounded to the C atom at the origin. The Hamiltonian of the impurity in the Floquet representation is and the hybridization term is Note that the the coupling matrix element V does not depend on the radiation field as we are considering normal incidence, hence the phase factor appearing in Eq. (35) is zero. The vacancy limit can be obtained from here by taking V → ∞. We define the Green function matrix G with elements given by G ij = ⟨⟨f i , f † j ⟩⟩. Using the Dyson equation it can be written as where G nm (ω) = ∑ k G nm (ω, k) and G nm (ω, k) = ⟨⟨a nk , a † mk ⟩⟩. Explicit expressions for the latter propagators are and with The propagator G 11 (ω, k) can be obtained from G 00 (ω, k) by the substitution ω ↔ (ω− ̵ hΩ) while G r 10 (ω, k) = G a 01 (ω, k) * where r and a denote retarded and advanced, respectively.
The energies of the bound states (if they exist) are determined by the poles of the trace of Eq. (39). This can be found numerically (as it is done below) but to grasp the main physical ingredients it is better to analyze the problem perturbatively. The imaginary part of the retarded self-energy V 2 G r 00 (ω) is proportional to the LDOS of the irradiated pristine graphene projected onto the m = 0 Floquet subspace and has a dynamical gap centered at ̵ hΩ 2. Its real part, on the other hand, is non zero inside the gap and diverges at the gap edges with different signs on each edge. As a consequence, to the lowest order in the impurity hybridization, the impurity spectral density (∝ −Im(G r 00 (ω))) has always a pole within the dynamical gap with an energy given by ω − − V 2 G r 00 (ω) = 0. Assuming, for the sake of argument, that = 0 , it is easy to see that in the same order and in the m = 1 Floquet subspace there is a bound state symmetrically positioned with respect to the gap center.
These results are in fact exact since G 01 (ω) = G 10 (ω) = 0 within the dynamical gap-we checked this numerically (see Fig. 10) but it can also be obtained from Eq. (41) in the low energy limit where φ k (A k , D(ω, k)) is odd (even) under the change k → −k. Therefore, there are two bound states, belonging to the m = 0 and m = 1 Floquet replicas, whose energies are given by the zeroes of ω − − V 2 G 00 (ω) and ω − − ̵ hΩ − V 2 G 11 (ω), respectively.  Figure 11 shows a color map of the local Floquet spectral density (corresponding to the three sites around the adatom) calculated using the Chebyshev method, described in Sec. V, as a function of the hybridization matrix element V for different values of . We found that while the energies of the bound states depend on the energy of the adatom, these states are always present regardless of the size of the hybridization. The symmetry between replicas is broken if ≠ 0 and it is only recovered in the limit of very large hybridization where the problem reduces to that of a vacancy. In this vacancy limit (V → ∞), the position of the bound states, are given by the solution of G r 00 (ω) = 0 and G r 11 (ω) = 0 (indicated by the arrows in Fig. 10), being the spectrum within the dynamical gap symmetric with respect to the gap center.
Interestingly, when looking at the weight of each of these states on the adatom and the three carbon atoms around it, one finds that they belong to a single replica. This particular result is a consequence that the coupling between the adatom and the layer of graphene was considered unaffected by the radiation field-see Figure 11.

VII. CONCLUSIONS
In summary, we have presented a detailed study of the Floquet bound states associated to defects in graphene illuminated by a laser. In particular, we focus on the bound states at the dynamical gap ( ̵ hΩ 2) using both analytical and numerical techniques applied to different defect types.
On one hand we consider large hole-like defects with different terminations. In this case, we show how the number of bound states increases with the defect radius and that the spectrum depends on the shape and type of lattice termination. In the case of cZZBC we proved analytically that in the limit of large radii the discrete bound states can be seen as nanoribbonlike chiral states 16,17 with a quantized linear quasi-momentum, as might have been anticipated. Staggered like potential (infinity mass boundary conditions) was also discussed with similar results, except that in this case there is a clear distinction between the two Dirac cones, and only one of them support chiral bound states. The chiral nature of the states was corroborated by an explicit calculation of the probability currents around the defect in the two analytical cases we presented.
On the other hand, we also consider point-like defects such as vacancies and adatoms and show that they also exhibit bound states around them. While the bound states spectrum depends on the value of the adatoms' orbital energy ( ) in the large hybridization or vacancy limit, they remain close to the bottom (top) border of the gap in the m = 0 (m = 1) replica.
Following the argument presented in Ref. [19] one can anticipate that additional bound states will also appear inside the high order gaps induced by high order photon processes. The contribution of such states to the spectral density projected onto the m = 0 replica is parametrically smaller provided η ≪ 1.
It remains a challenge for future work to evaluate the effect of these bound states on the bulk transport properties of dirty samples.

Appendix A: Boundary conditions
As we already mentioned in Sec. III, an arbitrary BC can be imposed by knowing the matrix M and their action on the wavefunction evaluated at the boundary: Ψ = M Ψ [58]. It can be demonstrated that boundary conditions are determined by two unit vectors:ν acting on the isospin (valleys) andn acting on the pseudospin (sublattices) [59]. In the isotropic representation M = (ν ⋅ τ ) ⊗ (n ⋅ σ), where τ and σ are the Pauli's matrices belonging to the isospin and pseudospin subspaces, respectively. In the following, we show the explicit form of the matrix M for regular polygons, included the circle as the limit case, and the three kinds of BCs considered in this work. For both, ZZBC and ABC/IMBC,n = ±ẑ (the sign depends on the sublattice termination) andn(ϕ) =ẑ ×n B (ϕ), respectively (see Fig. 12). In the latter expression,n B (ϕ) is the normal unit vector located at the edges of the defects pointing outward from the region of interest-for our purpose, this unit vector pointing to the center of the defects. For simplicity, we introduce the angle γ p related to the pseudospin degree of freedom. Thus, we can handle both types of boundary conditions at the same time by writinĝ and chose γ p = 0(π) or γ p = π 2 in order to select one or another type of BC. It must be noted that while z-component is exclusively related with the ZZBC, the xy-components are related with ABC and IMBC-the difference between two latter types of BCs resides in the isospin ν i.e., in the details of the lattice terminations. For a regular polygon with N sides, the normal unit vector pointing inwards has the form . It is useful to rewrite this quantity as a Fourier series where A m,N = sinc(π N ) (1 − mN ) and sinc(x) = sin x x.
It is straightforward to see that for circular defects we have lim N →∞ A m,N = δ m0 .
Analogously, for the isospin degree of freedom:ν =ẑ and ν ⋅ẑ = 0 for the ZZBC/IMBC and ABC, respectively. Introducing now the angle γ i , we can write this all three BCs in the formν where γ i = 0 for both, ZZBC and IMBC-for these BCs K and K ′ cones are decoupled. For the ABC however,ν lies on the xy plane, i.e., γ i = π 2-the Φ-phase is only relevant for the ABC, however, the analytic solutions of the ABC is out of the scope of this work. Finally, the matrix M in terms of the angles (γ i , γ p ) is and the analogous to the set of conditions (20), is The dependence of M with the polar angle ϕ relies on the pseudospin contribution. Triangles and hexagons are the unique regular polygons with well defined zigzag terminations. Therefore, the angle γ p for the ZZBC can behave in two different ways: it can be constant along the boundary of the defect (triangular defects), or it can alternate between 0 and π depending on the sublattice termination (hexagonal defects) (see Fig. 12). In order to tackle circular defects with ZZBCs, one is tempted to define the circle case as the limit of a polygon with a N large enough and an alternating n = ±z on their faces, corresponding to different sublattices terminations. However, this artificial limit is misleading because of is not possible construct such a defect, i.e., a regular polygon with N > 6 whose edges were constructed exclusively of zigzag neither armchair terminations. For simplicity, throughout this article we only work with ZZBC for triangular defects, in such a way that M is ϕ-independent. In this case, introduce the first condition of the set (A7) into Eq. (A6) leads to ψ B,l (ϕ, ξ 0 ) = ψ ′ B,l (ϕ, ξ 0 ) = 0 -we emphasize that, in the isotropic representation, ψ = (ψ A , ψ B , −ψ ′ B , ψ ′ A ) T must be used. Thus, there are two equations per cone [Eqs. (22) for the K cone], one for each Floquet replica, which allow us to find the relation between coefficients c + and c − , and then, the quasi-energies µ l [solutions of the Eqs . (23) and (25)].
On the other hand, for the ABC and the IMBC, the dependence of matrix M with the polar angle ϕ can not be avoided whatever the number of sides of the polygon considered. Even in the limit of circular defects: lim N →∞ Ξ N (ϕ) = ie −iϕunlike the ZZBC and the ABC, the circular defect is well defined for the IMBC because of this kind of BC does not depend on the details of the terminations at the edges (zigzag, armchair or mixing of them). As a consequence, for the IMBC the strategy to find the quasi-energies is quite different from that of the ZZBC (see App. C). For the K cone, the Floquet state restricted to n = 0 and n = 1 Floquet subspaces, has the form where the components are Hence, We also notice that λ − = λ * + and K ν (z * ) = K * ν (z). Because of cZZBC and IMBC do not mix different valleys [see Eq. (A5)], we can impose normalization conditions for each valley in an independent way. According to Eqs. (B1), and the angular dependence of the components of the Floquet state given by (B2), the normalization constant results timeindependent (1 + µ l ) 2 ξ dξ.
In order to obtain solutions belonging to the K ′ cone, the isotropic representation requires that: ψ A,l → −ψ ′ B,l and ψ B,l → ψ ′ A,l . By doing these replacements, the same procedure applied in Sec. II leads to a set of equations analogous to Eqs. (12)-and their respective boundary condition ψ ′ B (ξ 0 ) = 0which, in principle, must be solved again. However, for the cZZBC case, latter set of equations and their boundary conditions can be obtained from that of belonging to the K cone by doing follow changes: (µ, l) → (−µ, −l). Doing so, for the K ′ cone we have (B7) The time averaged probability density current (over one period) only has an angular component as it is shown in Sec. IV. Then, the density currents for both cones are J l = −2 Im{e iϕ u 1A,l u * 1B,l + u 0A,l u * 0B,l }, J ′ l = −2 Im{e −iϕ u ′ 1A,l u ′ * 1B,l + u ′ 0A,l u ′ * 0B,l }. (B9) Hence, it is straightforward to see that J l = J ′ −l . On the other hand, there is no any transformation between the K and K ′ cones for the IMBC case which simultaneously leaves invariant the set of differential equations and their respective boundary condition. Therefore, the set of quasienergies for the K ′ cone must be founded following the same procedure used for the K cone. The isotropic representation imposes that (B10) For the K ′ cone, coefficients c + and c − are related now by the phase e iθ = ω ′ + β + (ω ′ − β − ), with ω ′ ± = (1 + µ l )K l (λ ± ξ 0 ) − λ ± K l+1 (λ ± ξ 0 ). The time averaged probability density currents for each cone are also given by Eqs. (B8) and (B9). Nevertheless, there is no any relation between J l and J ′ l .
Appendix C: Solutions for the IMBC -Polygonal defects.
For simplicity, we will only tackle the IMBC, which does not mix cones. In this case, introducing the third condition of the set (A7) into Eq. (A6) leads to a mixing of solutions with different l quantum numbers due to the aforementioned dependence, i.e. where the upper (lower) sign refers to the K(K ′ ) cone and components given by Eq. (B2), [(B7)] were used. We also have to account for the dependence of the coordinates of the edges with the polar angle ϕ, i.e. ξ 0 (ϕ). For regular polygons with N sides, the points located at the edges can be written as where coefficients a m,N are given by R 0 is the apothem of the polygon andR 0 = R 0 a 0,N represents the mean value of their radii. For triangles and hexagons, a 0,3 = 3 ln(2 + √ 3) π and a 0,6 = 3 ln 3 π, respectively. In the large N limit, the deviations of R(ϕ) with respect toR 0 are small and we can expand the modified Bessel functions of second kind K ν appearing in Eqs. (C1) to first order on the deviation. That is, K ν (λξ 0 (ϕ)) ≃ K ν (λξ 0 ) + ∂K ν (λξ 0 ) ∂ξ 0 ξ 0 ξ 0 (ϕ) −ξ 0 . It is straightforward to see that only for circular defects, the mixing among different l quantum numbers is removed, since lim N →∞ a m≠0,N = lim N →∞ A m≠0,N = 0 and lim N →∞ a 0,N = lim N →∞ A 0,N = 1.
Finally, in order to find the quasi-energies, the infinite series in the Eqs. (C5) must be truncated. Doing so, it is possible to write a system with 2d equations for d quasi-energies (each quasi-energy introduce two additional coefficients: c + and c − ) and then, find their solutions.