{'fraction_non_alphanumeric': 0.028921782827179374, 'fraction_numerical': 0.010924186148131964, 'mean_word_length': 4.503832857357583, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear Response of Bilayer Graphene at Terahertz Frequencies

A density-matrix formalism within the length gauge is developed for the purpose of calculating the nonlinear response of intrinsic bilayer graphene at terahertz frequencies. Employing a tight-binding model, we find that interplay between the interband and intraband dynamics leads to strong harmonic generation at moderate field amplitudes. Specifically, we find that at low temperature (10 K), the reflected field of undoped suspended bilayer graphene exhibits a third harmonic amplitude that is 30% of the fundamental in the reflected field for an incident 1 THz single-cycle pulse with a field amplitude of 1.5 kV/cm. More interestingly, we find that as the central frequency of the incident radiation is increased, the third harmonic amplitude also increases; reaching a maximum of 53% for an incident frequency of 2 THz and amplitude of 2.5 kV/cm.


I. INTRODUCTION
Commencing in 2004 with the first observation of free standing graphene -along with the measurement of its electronic properties [1]-there has been an abundance of both experimental and theoretical work performed on graphene and its multi-layered counterparts [2][3][4][5][6][7]. The impetus behind this research arises from the intriguing electrical, mechanical, and optical properties possessed by graphene. Many of the intriguing properties of monolayer graphene (MLG) are also present in bilayer graphene (BLG). These include a zero-gap energy dispersion, a very high carrier mobility, high thermal conductivity, and the ability to tune electric properties by adjusting the carrier density through gating [2,8].
Despite their similarities, BLG differs from MLG in a number of ways. One important distinction is found in the density of states of both systems. Close to the Dirac points -the points within the first Brillouin zone at which the conduction and valence bands touch -one can show that although the interband matrix elements for both systems are essentially the same, in BLG, the number of states within 10 meV of the Dirac point is approximately 30 times larger than in MLG. This difference in the density of states will have an effect on the interband transitions of the two systems. In order to probe this energy range one needs to probe BLG with radiation that has frequencies in the terahertz range.
As with MLG, the absorption of optical and THz radiation by BLG can be understood in terms of interband and intraband transitions. The energies of photons at THz frequencies are ideally suited to the study of intrinsic (undoped) BLG, where the Fermi level is at the Dirac point. At these frequencies, interband transitions can take place near the Dirac point, while the injected electrons and holes are subsequently strongly driven within their respective bands. As we shall see, the ability of THz radiation to induce both intraband and interband dynamics is fundamental in the emergence of a strong nonlinear response in BLG.
Earlier studies on MLG have suggested the presence of a strong nonlinearity at both optical and THz frequencies [9][10][11][12][13]. Theoretical studies of BLG have indicated a strong nonlinear effect in the THz to far-infrared regime, whereby a moderate electric field can result in third harmonic generation at room temperature [14]. Experiments have been performed with the intent of observing the nonlinear THz response of graphene, specifically third harmonic generation. Although it has been observed by using a 45-layer sample, it has not yet been successfully observed in MLG or BLG [6,7]. Almost all experimental investigations of the nonlinear response of graphene have observed a response that is dominated by intraband dynamics, due to the high doping level of the samples. The large Fermi energy of the doped system diminishes the interband current due to Pauli blocking. For more moderate Fermi energies -tens of meV -the interband dynamics can make a contribution that is dependent on the field amplitude. In recent theoretical work on MLG, it was shown that if the Fermi level is reduced to within a few meV of the Dirac point, the magnitude of the interband current becomes comparable to the intraband current, resulting in the presence of a strong nonlinearity [15]. Specifically, it was found that odd harmonics of the THz field should be generated. As this interplay between interband and intraband dynamics is a result of the zero-gap band structure, we also expect to see a similar interplay in the carrier dynamics in BLG.
In this paper, we employ a four-band tight-binding method to model the intraband and interband dynamics of undoped suspended bilayer graphene in response to a single-cycle THz pulse in the range of 1 to 5 THz. We use this model to calculate the dependency of the third harmonic response of BLG on the graphene temperature and the central frequency of the pulse. The current densities and the generated harmonics are calculated numerically for two values of the temperature, 10 and 100 K. We find that the ratio of the amplitude of the third harmonic to the fundamental in the reflected field is reduced by an order of magnitude as the temperature is increased from 10 to 100 K. Finally, we examine the nonlinear response at a number of central THz frequencies. We find that as the central frequency increases from 1 to 5 THz, there is an increase in the ratio of the third harmonic to the fundamental; reaching a maximum at a central frequency of 2 THz.
The paper is organized as follows. In section II, we present the derivation of our theoretical model. A tightbinding method is used to obtain dipole matrix elements within the length gauge. The matrix elements are then used in a four-band density matrix formalism in order to calculate the interband and intraband current densities. In section III, we present both the linear and nonlinear THz response of numerical simulations for undoped BLG at different temperatures and central THz frequencies.
The conclusions are presented in section IV.

II. THEORY
The calculations that we perform are based on a theoretical approach employing a density-matrix formalism in the length gauge. A nearest-neighbor tight-binding model is used to treat the π-electrons in the graphene, which are taken to provide the conduction electrons only [16].

II.1. Energy Bands
The tight-binding expression for the Bloch states is given by where A n (k) is a normalization factor, n labels the conduction and valence bands, and the sum is over the Bravais lattice vectors R. The sublattice coefficients, C n i (k), are associated with the four carbon atoms within the unit cell; the ϕ pz (r) are the 2p z orbitals of carbon. The index i indicates a sum over the basis vectors r A1 , r B1 , r A2 , r B2 , which give the position of sublattice sites A 1 and B 1 in the top layer, and A 2 and B 2 in the bottom layer (see Fig. 1). Explicitly, they are given by The vector of the sublattice coefficients is given by Within the basis of these vectors, the secular equation for BLG can be expressed as with the eigenvalue of band n given by ǫ nk .
The intralayer hopping energy, t , and the interlayer hopping energy, t ⊥ , are approximately equal to 3.03 eV and 0.3 eV, respectively [17]. The function f (k) ≡ (1 + e −ik·a1 + e −ik·a2 ), is a result of the nearest neighbor intralayer electron hopping, where the a i are the primitive translation vectors of graphene, given explicitly by Here a o is the nearest-neighbour separation (a o ≃ 1.42Å), not the length of the primitive vectors. Equation (2) can be solved exactly to yield expressions for the eigenvalues of BLG. In order of increasing energy, these are given by where the first two expressions represent the valence bands, the second two the conduction bands and where we defineǫ (k) = t 2 ⊥ + 4|f (k) | 2 . Due to the symmetry between the sublattices, the conduction and valence states in graphene are degenerate at two Dirac points, given by: For energies within a few hundred meV of the Dirac points, we can expand the crystal momentum around the Dirac points as k = K + δk and k ′ = K ′ + δk, where δk = k xx + k yŷ . With this expansion, we can express the band energy near both Dirac points asǫ (k) as where k is the magnitude of the crystal momentum k ≡ |δk|, α = 4 2 v 2 F , where v F = 3a o t ⊥ /2 is the Fermi velocity. The energy dispersion given in Eq. (4) for the four bands is shown in Fig. 2 for k-vectors near the K-Dirac point. Using the expressions for the eigenvalues in Eq. (4), we may also solve for the corresponding eigenstates exactly, to obtain where e iχ(k) = f (k) / |f (k)|, and we have supressed the explicit k-dependencies for simplicity.
As discussed in previous studies [15,18,19], modelling the nonlinear response of semiconductors in a limited band model is only expected to be accurate if one uses the length gauge. If instead the velocity gauge is employed, unphysical low-frequency divergences arise in the nonlinear response. As these divergences affect the response in the THz regime, one must develop sum rules in order to remove them. Moreover, trying to obtain the response to higher and higher order in the field, these sum rules become analytically intractable. Consequently, when using a limited basis of bands, the length gauge is advantageous. Hence, to model the THz interaction of BLG, we employ the length gauge Hamiltonian, expressed as H = H 0 − er · E, where H 0 is the full Hamiltonian of unperturbed BLG, e = −|e| is the charge of an electron, r is the electron position vector, and E(t) is the THz electric field at the graphene. Given that we only consider normally incident plane waves, the field is uniform over the graphene sheets.
The carrier dynamics in BLG are calculated by solving the equations of motion for the density matrix in the basis of the Bloch states describing the four bands of BLG. For these, we require the matrix elements of the Hamiltonian between the various Bloch states: The matrix elements of r between Bloch states can be shown to be given by [15,19] n, k |r |m, where the connection elements are defined as where Ω is the volume of a unit cell and u nk (r) is the periodic part of the Bloch function. These connection elements have been calculated using our tight-binding wave function. Ignoring the overlap of atomic wave functions on different atoms, near the K-Dirac point the interband connection elements can be shown to be given by where δk = cos(θ)x+sin(θ)ŷ andθ = −sin(θ)x+cos(θ)ŷ are, respectively, the radial and angular unit vectors in cylindrical coordinates with the origin at the K-Dirac point. Around the K ′ -Dirac point, we find that all the connection elements change sign except for those with components in thek-direction. Explicitly, we have Similarly, we can show that all the intraband connection elements -which are identical to the Berry connections of the respective bands [20,21] -are zero, i.e.
This is in agreement with the idea that the Berry phase of BLG does not have to be simply 2π, but rather can be an integer multiple of 2π [22]. The vanishing Berry connection leads directly to a zero Berry curvature for each of the four bands; which is expected due to inversion symmetry present in BLG. The dynamic equations for the density matrix elements can be found using the approach of Aversa and Sipe [19], which was used in similar work on MLG [15]. These equations are given by )/ , and ρ eq nm (k) = f n (k, T )δ nm is the carrier population in equilibrium when n = m, and is zero otherwise; where f n (k, T ) is the Fermi-Dirac distribution at temperature T . In our numerical work, we model the vacancy populations rather than the valence band electrons to allow us to only include states near the Dirac point, which greatly reduces computation time.
Previous studies have shown that the scattering times in graphene are on the order of tens of femtoseconds [6,7], therefore to accurately model the THz response we must take into account scattering processes. At the low carrier densities considered in this work, carrier-carrier scattering is expected to be relatively unimportant, and the dominant scattering processes will be defect scattering and electron-phonon scattering. To account for these mechanisms in our model, we treat scattering phenomenologically. We introduce an interband decoherence time, τ nm , for the interband coherences, ρ nm (k), where n = m. We assume the decoherence time to be independent of k. The populations, ρ nn (k), we take to relax back to Fermi-Dirac thermal distributions, f n (k, T ), with relaxation times, τ n . Since the interaction with THz radiation induces interband transitions, as the simulation proceeds, we adjust the temperature of the Fermi-Dirac distribution so that the carriers relax to the instantaneous carrier populations. Also, we neglect interband relaxation as it has been found that the intraband scattering times are much shorter than the time taken for conduction band electrons to relax to the valence band [23].
We employ a direct computational approach to solve the above equations, wherein we put k on a grid and step through time using a Runge-Kutta algorithm. To facilitate this, we make use of balanced difference approximations to the gradients. Now that we have the dynamic equations for the density matrix elements, we turn to determining the expression for the current density. Following the formalism of Aversa and Sipe [19] the current density can be expressed as Using the fact that p m = 1 i [r, H] , and decomposing the position operator into intraband and interband parts, r = r i + r e , we can write this as [19] where the trace is over the single electron states, and ρ(t) is the reduced density matrix with matrix elements ρ nm (k). The decompostion of the position operator allows us to define the total current density as the sum of an intraband contribution, J i , and an interband contribution, J e . Using our expression for the Hamiltonian, as well as the matrix elements of the position operator (Eq. 12), one may determine expressions for these contributions. The procedure is similar to that presented in recent work on MLG [24]. After considerable work, the intraband current density near the K-Dirac point can be shown to be given by while the interband current density is given by The sums over k in the current density expressions are restricted to a region near the K-Dirac point, and we include a factor of 2 to account for spin degeneracy. We also need to take into account the current density near the K ′ point. Similar to MLG, we find that due to the symmetry of the Brillouin zone, the current densities around both Dirac points are identical. To obtain the total current density of BLG, we then only need to multiply the results calculated at the K-point by two in order to account for this degeneracy.
In what follows, we consider a suspended BLG sample, and use the time-dependent current densities to calculate the transmitted and the reflected THz fields, using a procedure identical to that used for MLG [24]. We show in the next section that for low amplitude fields, we obtain the expected linear conductivities, which relate the THz field at the graphene sheet to the induced current densities. Furthermore, we have verified convergence in the nonlinear regime by altering the grid density, the extent of the grid, the time-step tolerance, and the polarization of the incident field.

II.2. Linear Response
In order to calculate the linear response of BLG to an incident field, we need to calculate expressions for the density matrix elements to first order in the THz field. Once we have these, we may then use Eq. (20) to express the first order intraband current density as whereǫ ′ is a shifted energy, given byǫ ′ =ǫ− t ⊥ . C(ω p ) is a frequency-dependent coefficient that includes the electric field: where we take our field to be harmonic, In Eq. (22), the sum is over the Fermi-Dirac distributions of the electrons and holes in the four bands of BLG. This is given expilicity by where β = k B T and µ is the chemical potential of the system.
As a check, we may evaluate the first order intraband current density in the limit of zero interlayer coupling. In this limit, t ⊥ → 0, only the integral of the sum over the distributions remains. Thus our intraband current density reduces to From here it is easy to read off the first order intraband conductivity as being This is the previously found result for the monolayer reponse [10], multiplied by two to account for the two layers of BLG.
For BLG (t ⊥ = 0), the integral in Eq. (22) must be evaluated numerically in general, however in the limit that T → 0, it may be evaluated analytically. In this limit the Fermi-Dirac distributions will behave as Heaviside step functions; having the effect that only energies lower than the Fermi level, E F , contribute to the integral. Assuming that our Fermi level is non-zero and located in the conduction bands, we only get contributions from the ρ c1c1 (k) and ρ c2c2 (k) distributions. In this case, we obtain two distinct contributions to the intraband current density. Thus, we are able to determine the intraband conductivities of both conduction bands in the limit of zero temperature. For the low-energy band c 1 , this is given by while for the high-energy band c 2 , we have Here the step function assures that the Fermi level must be greater than t ⊥ for ρ c2c2 (k) to provide a contribution.
Similarly, we may now use Eq. (21), and the connection elements from Eq. (14), to express the first order interband current density as where D(ω p ) is a frequency dependent coefficient given by and the function R(x) is dependent on the temperature, Fermi level and energy of the system, which is given explcitly by .
We may evaluate the above integral analytically in the limit of zero temperature and scattering, where T → 0 and τ → ∞, respectively. In this limit, the complex factors in each term reduce to a Dirac delta function -when we neglect the contribution of the principal part of the integral. Taking only the positive frequency portion of the above expression, we obtain for the first order interband conductivity (29) This expression agrees with those found in the literature [25,26]. We see that BLG exhibits a minimum conductivity that is associated with interband transitions. For very large incident frequencies, this minimum conductivity is given by |e| 2 2 , which is twice that found in monolayer graphene, as expected.
We may also use Eqs. (22) and (28) to determine the first order current densities at non-zero temperatures and scattering times, which is a closer approximation of experimental conditions. We present the results of this analysis in the next section.

III. SIMULATION RESULTS
We have shown that our analytic expressions for the linear THz response agrees with the literature in certain limits. Going beyond the linear regime to determine the response at higher order in the THz field is our main interest. We want to be able to examine the dependency of the THz response on the temperature and incident frequency of BLG. To achieve this aim we employ a simulation, which offers the ability to control the desired parameters, and calculate the interband and intraband current densities at the desired field amplitudes. The field transmitted from the BLG is calculated as a function of the current densities and the incident field. The transmitted field is then spectrally analyzed to determine its frequency components; high harmonic generation in the spectral composition indicating nonlinear behaviour.
We begin with a check on the first order calculations by comparing the linear conductivity derived from the first order equations of the previous section to that calculated by the full simulation. We then proceed to examine the higher order response of BLG, determining at which incident field amplitudes we might expect to see the largest generation of the third harmonic.

III.1. Linear Results
As an initial check, we compare on a single plot the conductivity arising from the full current density (intraband and interband), as calculated by both the simulation and the closed form expressions (Eqs. (22) and (28)). We plot these results in Figs. 3 and 4 for Fermi levels of 60 and 360 meV respectively, which correspond to energies of 0.2t ⊥ and 1.2t ⊥ . These values also allow us to make comparisons to similar work found in the literature [26]. In both instances, the scattering time is 80 f s, and the temperature is 50 K. As can be seen the agreement between the simulated and first order results is excellent.
In both plots we see features at incident freqeuncies of f = 0 and f ≃ 75 THz. The zero frequency peak is a result of the intraband Drude response. This feature vanishes for zero Fermi level and temperature, as the intraband transitions become negligible in this limit. The second peak at f ≃ 75 THz is due to the difference in energy between the two conduction bands of BLG (≃ 300 meV ). Specifically, it arises from the perfect nesting between the dispersions of the higher and lower conduction bands, c 2 & c 1 . These two bands sit directly on top of one another and differ in energy only by the constant displacement t ⊥ . Thus, the nesting of the bands leads to a large spectral weight at frequencies close to t ⊥ , which is approximately equivalent to 75 THz. The additional features present in the plots reflect other possible interband transitions, such as v 1 → c 1 or v 1 → c 2 . These transitions require incident radiation of energy t ⊥ and 2t ⊥ , respectively, and their spectral weights are largely dependent on the Fermi level.

III.2. Nonlinear Results
We now present our simulation results for the nonlinear THz response of BLG. In each of our simulations for BLG we took the sample to be undoped (µ = 0). The response was then calculated at field amplitudes ranging from 0.5 kV/cm to 2.0 kV/cm at two different temperatures: 10 K and 100 K. This is followed by the calculation of the response at several incident frequencies at a temperature of 10 K. For the system and fields modelled in this section, only the v 2 and c 1 bands contribute to the THz response.
Our input THz field is a sinusoidal Gaussian pulse with central frequency of 1 THz and full width at half maximum (FWHM) of 1 ps. Mathematically, this may be expressed as where t 0 and T F W HM are the temporal shift and full width at half maximum of the Gaussian pulse, respectively. The central frequency of the pulse is given by f 0 . We take our scattering and relaxation times to be fixed at 50 fs, which is an average of theoretical and measured values for such (sample-dependent) constants. This is a conservative estimate for undoped BLG at low temperature; an increase in this time will lead to stronger harmonic generation [15]. In the simulation, we keep the Fermi level at the Dirac point and as carriers are injected we raise the temperature of the distribution to which the carriers relax to account for the increase in the carrier density.

T=10 K Results
At a temperature of 10 K, the intrinsic thermal carrier density in BLG is approximately 1.4 × 10 10 /cm 2 and the electrons have an average energy of 1.03 meV. In comparison, at this temperature MLG has an intrinsic carrier density of approximately 2.0 × 10 8 /cm 2 , and an average electron energy of 1.89 meV. Thus, we find a much larger thermal carrier density in BLG than we do in MLG. However, the average energy of the carriers is slightly larger in MLG. In both cases, the average electron energy is much less than the central photon energy of 4.14 meV associated with a 1 THz pulse, thus to first order in the field, interband transitions are largely unaffected by the thermal energy.
All of the carriers (injected and thermal) are driven within the conduction and valence bands by the applied electric field. The subsequent motion leads to a blocking and unblocking of the electronic states that are available for extra carrier injection. This results in an interplay between intraband and interband dynamics that is paramount in producing the nonlinear response [15].
We now present results for a number of different field amplitudes. In Fig. 5, we plot the intraband and interband current densities for four different incident field amplitudes (0.5, 1, 1.5 and 2 kV/cm). All current densities are normalized to the peak value, E 0 , of the incident field such that, if the response were linear, these relative currents would be unchanged by an increase in incident field. This allows for a comparison of the currents at each field amplitude, and for the clear identification of any nonlinear behaviour. In what follows, we refer to these as relative current densities. The calculated relative intraband and interband current densities at these specific field amplitudes are shown in Figs. 5(a) and 5(b), respectively.
Looking first at the intraband current density, we find that as the field amplitude is increased, the relative intraband current decreases. This is due to the 'clipping' phenomena that is a property of linear dispersions, such as that found in MLG. Although BLG has a parabolic dispersion near the Dirac point, as we move away from this point, the dispersion becomes linear (see Fig. 2). Thus, the high field amplitudes drive carriers far enough away from the Dirac point such that they occupy the linear part of the band structure. In this region, the carrier velocities become constant and the field is not able to drive them to higher velocities; the intraband current is clipped as a result.
We note, however, that the current clipping nonlinearity in the intraband current density is much smaller than the nonlinearity seen in the interband current. We also note that due to the large intrinsic carrier density of BLG due to the large number of states, carrier injection has only a minor effect on the intraband current density at these field amplitudes. This is in contrast to MLG, where one finds a substantial increase in the relative intraband current density arising from the injection of carriers [15].
Examining next the interband current density, we observe large distortions for all of the field amplitudes. Based on the strength of these distortions, we expect the interband current to be a large source of nonlinearity at these field amplitudes. The motion of the charge carriers in the conduction band gives rise to an interplay between intraband and interband dynamics. We can visualize this interplay between interband and intraband motion via a plot of the electron density in the conduction band, as shown in Fig. 6. As the pulse hits the BLG sample, the population of carriers in the conduction band is driven away from the Dirac point by the incident THz field; this driving of the carriers within each band is the source of the intraband current. Fig. 6(a) shows the initial electron density before the pulse hits. In Fig. 6(b) we see the density at t = 1.75 ps, where the magnitude of the amplitude of the incident field is a maximum. We can see the distribution has been driven to the right of the Dirac point. This intraband motion also opens up electron states near the Dirac point that were previously occupied before the pulse arrived. As a result of the reduction in Pauli blocking, carriers from the valence band v 2 may be injected into the conduction band c 1 near the Dirac point. This is evidenced in the dark, high-density regions that appear above and below the k y = 0 line to the right of the Dirac point in Fig. 6(b). This carrier injection is the source of the interband current.
Thus, we see that there is an interesting relationship between the intraband and interband dynamics; the motion of carriers within each band has an effect on both interband and intraband current densities. Importantly, this relationship manifests itself in the nonlinear response of BLG.
We now plot the time-dependent reflected fields normalized to the peak value of the different incident field amplitudes in Fig. 7(a), and the spectral responses normalized to the peak amplitude at the fundamental frequency (1 THz) in Fig. 7(b). In the upper plot we can see strong distortion of the reflected fields for all of the selected field amplitudes. Note that because there is no substrate the induced current is the entire source of the reflected field, which only has a peak amplitude of approximately 2 V/cm for the incident field of 1 kV/cm. To see harmonic content in the reflected field, Fig. 7(b) we plot the Fourier transform of the reflected fields normalized to the peak reflected field at the fundamental frequency. We find that at the lowest field amplitude of 0.5 kV/cm, there is a clear third harmonic signal. Furthermore, we see a maximum in the third harmonic of the relfected field at an incident field amplitude of about 1.5 kV/cm. At this amplitude the third harmonic reaches a maximum of approximately 30% of the reflected spectral peak at the fundamental. As we increase the incident field amplitude further, we begin to see a decrease in the third harmonic.
Moreover, at these high field amplitudes we also see the presence of the 5th harmonic, which is largely absent at lower field amplitudes. Its amplitude takes a maximum value of approximately 15% the spectral peak at the fundamental for an incident field of 2 kV/cm. The appearance of the fifth harmonic indicates that at the higher field amplitudes fifth order processes are significant. Just as was found with MLG [15], the fifth order response results in a decrease in the third harmonic.

T=100 K Results
We next present results for simulations at a temperature of 100 K. At this temperature, the intrinsic thermal carrier density is 3.4 × 10 10 /cm 2 and the average thermal electron energy is 10.8 meV. This is compared to the average photon energy of 4.14 meV associated with a 1 THz pulse. Thus, the Pauli blocking of the interband transitions is increased due to the increased thermal populations of carriers; reducing the overall interband current density. Additionally, the increase in thermal carriers gives rise to an increase in the intraband current density. Thus, at this temperature the intraband current is dominant, resulting in a reduction in the interplay between intraband and interband dynamics, and ultimately a reduction in high harmonic generation.
We can see this by looking again at both the reflected field and the spectral response for each of the incident field amplitudes. The normalized time-dependent reflected fields for the different field amplitudes are shown in Fig. 8(a), and the spectral responses normalized to the peak amplitude at the fundamental frequency (1 THz) are presented in Fig. 8(b).
In the upper plot we see far less distortion in the reflected fields than we do at T=10 K, for all of the selected field amplitudes. The reduction in distortion is again due to the increase in the thermal carrier density and subsequent increase in Pauli blocking. In turn this results in diminished harmonic generation, which is shown in Fig. 8(b). We find a maximum in the third harmonic of the relfected field at an incident field amplitude that is now approximately 1 kV/cm. At this amplitude it reaches a maximum of only approximately 2.8% of the reflected peak at the fundamental, which is still greater than that found in MLG at the same temperature [24]. At higher field amplitudes, we also see the presence of the 5th harmonic. Its amplitude takes a maximum value of approximately 1.5% of the peak at the fundamental for an incident field of 2 kV/cm.
Thus, we find that as the temperature of the system is increased, the nonlinear response -specifically, high harmonic generation -is greatly diminished as a result of the reduction in the interplay between interband and intraband dynamics. Therefore, if one is to observe high harmonic generation (HHG) experimentally in BLG, it is apparent that low temperatures are necessary. This allows for the maximization of the dynamic interplay that is paramount to the nonlinear repsonse.

Effect of Central Frequency
Finally, we study the effect that the central frequency of the incident pulse has on the generation of the third harmonic. The interplay between intraband and interband dynamics has a strong dependence on the central frequency, because when the central frequency is increased, the carriers are injected farther from the Dirac point. In Fig. 9 , we plot the normalized third harmonic amplitude as a function of the incident THz field amplitude for central frequencies of 1.0, 2.0 and 5.0 THz. In each case, we adjust the duration of the pulse such that the product of the central frequency and the pulse duration is constant, so that the pulse remains single-cycle for each of the central frequencies considered.
As the central frequency of the incident pulse is increased, we observe a maximum third harmonic amplitude that is larger in comparison to the simulations performed at 1 THz. For a central frequency of 2 THz, the normalized third harmonic amplitude peaks at approximately 53% of the fundamental for a field amplitude of 2.5 kV/cm. This is compared to a value of 30% for the central frequency of 1 THz at a field amplitude of 1.5 kV/cm.
As the pulse frequency is increased further to 5 THz, the maximum normalized third harmonic is found to decrease to 40% of the fundamental. Furthermore, this maximum value is not obtained until the input field amplitude reaches 5 kV/cm.
The observed increase in the nonlinear interplay with increasing central frequency is a result of a number of factors. First, due to the approximately inverse relationship between the linear conductivity and the frequency, the intraband current density decreases as we increase the central frequency of the pulse. This in turn increases the ratio of the third harmonic to the fundamental, due to the reduction in the part of the reflected fundamental that arises due to the intraband current.
Additionally, we find that as the incident frequency is increased there is an increase in the interband current density. The increase arises because the interband transitions that are resonant with the 2 THz and 5 Thz photons are mostly unoccupied when the pulse arrives, whereas at lower frequencies they may be occupied by either thermal or injected carriers. This larger interband current density results in an increased dynamic interplay, and thereby yields a larger third harmonic field.

IV. SUMMARY
We have presented a detailed study of the nonlinear response of unbiased bilayer graphene at THz frequencies. A theoretical model has been developed, which is based on the dynamic equations of density matrix elements, employing the length gauge. The model enabled the calculation of eigenvalues and eigenvectors of BLG, as well as interband and intraband connection elements. Expressions for interband and intraband current densities were also derived.
Through the use of simulation, we determined solutions of the density matrix dynamic equations. These solutions were then applied to the study of high harmonic generation in BLG. Investigating the effect of system parameters -such as the central THz frequency and the ambient temperature -on third harmonic generation was the main focus.
Our results show that for a temperature of 10 K, scattering time of 50 fs and incident field of 1 THz, the third harmonic can be as large as 30% of the fundamental for an incident THz field of 1.5 kV/cm. As we increase the temperature to 100 K, we find that the maximum third harmonic generation is reduced by an order of magnitude, due to the reduction in interband transitions.
Finally, we found that as the central THz frequency is increased from 1 to 5 THz, we see an increase in the third harmonic amplitude; reaching a maximum of 53% of the fundamental for a 2 THz incident field at 2.5 kV/cm. These results may be of use to experimentalists aiming to probe the nonlinear response of bilayer graphene.
To experimentally observe the high harmonics we predict for BLG, one must consider the dynamic range of the THz spectrometers -defined as the ratio of the frequency dependent signal strength to the detected noise floor [27]. For a 1 THz incident field of 1 kV/cm, we find the peak amplitude of the reflected field to be approximately 2.2 V/cm (53 dB less than the incident field). Thus, a detection technique that allows for a dynamic range larger than 53 dB is required for the measurement of the reflected signal. Such a dynamic range is experimentally feasible; a dynamic range of 90 dB has been reported in recent work [28].