Wave packet dynamics in slowly modulated photonic graphene

TL;DR

This paper mathematically analyzes wave packet dynamics in photonic graphene with honeycomb structures, showing that wave packets near Dirac points evolve according to a Dirac equation with a varying mass, illuminating topological phenomena.

Contribution

It provides a rigorous mathematical derivation of wave packet evolution in modulated photonic graphene, linking electromagnetic wave behavior to Dirac equations with variable mass.

Findings

Wave packets near Dirac points follow a Dirac equation with varying mass.

The analysis offers mathematical insights into topological phenomena in photonic graphene.

Wave packet dynamics are governed by the spectral concentration at Dirac points.

Abstract

Mathematical analysis on electromagnetic waves in photonic graphene, a photonic topological material which has a honeycomb structure, is one of the most important current research topics. By modulating the honeycomb structure, numerous topological phenomena have been observed recently. The electromagnetic waves in such a media are generally described by the 2-dimensional wave equation. It has been shown that the corresponding elliptic operator with a honeycomb material weight has Dirac points in its dispersion surfaces. In this paper, we study the time evolution of the wave packets spectrally concentrated at such Dirac points in a modulated honeycomb material weight. We prove that such wave packet dynamics is governed by the Dirac equation with a varying mass in a large but finite time. Our analysis provides mathematical insights to those topological phenomena in photonic graphene.

Equations

$$\partial_{t}^{2}\psi-\nabla\cdot\big{(}W_{\varepsilon}(\boldsymbol{x})\nabla\psi\big{)}=0,\quad(\boldsymbol{x}\in\mathbb{R}^{2},~{}t>0)$$

$$W_{\epsilon }(\mathbf{x})=A(\mathbf{x})+\epsilon \kappa (\epsilon \mathbf{x})B(\mathbf{x}).$$

$${\mathcal{L}}^A{\Phi }_j(\mathbf{x})=E_D{\Phi }_j(\mathbf{x}),j=1,2.$$

$$\left\{\begin{split}\psi(\boldsymbol{x},0)=&~{}\varepsilon\Big{[}a_{1}(\varepsilon\boldsymbol{x})\Phi_{1}(\boldsymbol{x})+a_{2}(\varepsilon\boldsymbol{x})\Phi_{2}(\boldsymbol{x})\Big{]},\\ \partial_{t}\psi(\boldsymbol{x},0)=&~{}\varepsilon\Big{[}b_{1}(\varepsilon\boldsymbol{x})\Phi_{1}(\boldsymbol{x})+b_{2}(\varepsilon\boldsymbol{x})\Phi_{2}(\boldsymbol{x})\Big{]},\end{split}\right.$$

$$\left\{\begin{aligned} \psi(\boldsymbol{x},0)=&~{}\varepsilon\Big{[}\alpha_{10}(\varepsilon\boldsymbol{x})\Phi_{1}(\boldsymbol{x})+\alpha_{20}(\varepsilon\boldsymbol{x})\Phi_{2}(\boldsymbol{x})\Big{]},\\ \partial_{t}\psi(\boldsymbol{x},0)=&~{}i\sqrt{E_{D}}\varepsilon\Big{[}\alpha_{10}(\varepsilon\boldsymbol{x})\Phi_{1}(\boldsymbol{x})+\alpha_{20}(\varepsilon\boldsymbol{x})\Phi_{2}(\boldsymbol{x})\Big{]},\end{aligned}\right.$$

$$\displaystyle\psi(\boldsymbol{x},t)=e^{i\sqrt{E_{D}}t}\varepsilon\Big{[}\alpha_{1}(\varepsilon\boldsymbol{x},\varepsilon t)\Phi_{1}(\boldsymbol{x})+\alpha_{2}(\varepsilon\boldsymbol{x},\varepsilon t)\Phi_{2}(\boldsymbol{x})\Big{]}+\eta(\boldsymbol{x},t),$$

$$\left\{\begin{aligned} i\partial_{T}\alpha_{1}-\frac{v_{{}_{F}}}{2\sqrt{E_{D}}}\big{(}i\partial_{X_{1}}-\partial_{X_{2}}\big{)}\alpha_{2}+\frac{\vartheta_{\sharp}\kappa(\boldsymbol{X})}{2\sqrt{E_{D}}}\alpha_{1}=0\\ i\partial_{T}\alpha_{2}-\frac{v_{{}_{F}}}{2\sqrt{E_{D}}}\big{(}i\partial_{X_{1}}+\partial_{X_{2}}\big{)}\alpha_{1}-\frac{\vartheta_{\sharp}\kappa(\boldsymbol{X})}{2\sqrt{E_{D}}}\alpha_{2}=0\end{aligned}~{},\right.$$

$${\alpha }_1(\mathbf{X},0)={\alpha }_{10}(\mathbf{X}),{\alpha }_2(\mathbf{X},0)={\alpha }_{20}(\mathbf{X}).$$

$$\underset{0\le t\le \rho {\epsilon }^{{}^{-1}}}{\sup }\Vert \eta (\mathbf{x},t){\Vert }_{H^s({\mathbb{R}}^2)}\le C{\epsilon }^{1-\nu },$$

$$\mathcal{L}^{A}=-\nabla\cdot\big{(}A(\boldsymbol{x})\nabla\big{)},~{}\mathcal{L}^{B}=-\nabla\cdot\big{(}B(\boldsymbol{x})\nabla\big{)},~{}\mathcal{L}^{B}_{\kappa}=-\nabla\cdot\big{(}\kappa(\varepsilon\boldsymbol{x})B(\boldsymbol{x})\nabla\big{)}.$$

$$\chi (|\mathbf{k}|\le C)f(\mathbf{k}):=\{\begin{array}{r}f(\mathbf{k}),\text{if}|\mathbf{k}|\le C,\\ 0,\text{otherwise}.\end{array}$$

$${\sigma }_1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),{\sigma }_2=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right),{\sigma }_3=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right).$$

$$\big{\langle}f_{1}(\boldsymbol{x}),f_{2}(\boldsymbol{x})\big{\rangle}=\int_{\mathbb{R}^{2}}\overline{f_{1}(\boldsymbol{x})}f_{2}(\boldsymbol{x})d\boldsymbol{x},\qquad\big{\langle}g_{1}(\boldsymbol{x}),g_{2}(\boldsymbol{x})\big{\rangle}_{L^{2}(\Omega)}=\int_{\Omega}\overline{g_{1}(\boldsymbol{x})}g_{2}(\boldsymbol{x})d\boldsymbol{x}.$$

$$\hat{f}(\boldsymbol{\xi})=\int_{\mathbb{R}^{2}}e^{-i\boldsymbol{\xi}\cdot\boldsymbol{x}}f(\boldsymbol{x})d\boldsymbol{x},\qquad\tilde{f}_{b}(\mathbf{k})=\big{\langle}\Phi_{b}(\cdot;\mathbf{k}),~{}f(\cdot)\big{\rangle}~{}.$$

$$\Lambda =\{\mathbf{mv}=m_1{\mathbf{v}}_1+m_2{\mathbf{v}}_2:\mathbf{m}=(m_1,m_2)\in {\mathbb{Z}}^2\}=\mathbb{Z}{\mathbf{v}}_1\oplus \mathbb{Z}{\mathbf{v}}_2,$$

$$\Omega =\{{\theta }_1{\mathbf{v}}_1+{\theta }_2{\mathbf{v}}_2:0\le {\theta }_j\le 1,j=1,2\},$$

$${\mathbf{v}}_1=\left(\begin{array}{c}\frac{\sqrt{3}}{2}\\ \frac{1}{2}\end{array}\right),{\mathbf{v}}_2=\left(\begin{array}{c}\frac{\sqrt{3}}{2}\\ -\frac{1}{2}\end{array}\right).$$

$${\Lambda }^{\ast }=\{\mathbf{mk}=m_1{\mathbf{k}}_1+m_2{\mathbf{k}}_2:\mathbf{m}=(m_1,m_2)\in {\mathbb{Z}}^2\}=\mathbb{Z}{\mathbf{k}}_1\oplus \mathbb{Z}{\mathbf{k}}_2,$$

$${\mathbf{k}}_1=\frac{4\sqrt{3}}{3}\pi \left(\begin{array}{c}\frac{1}{2}\\ \frac{\sqrt{3}}{2}\end{array}\right),{\mathbf{k}}_2=\frac{4\sqrt{3}}{3}\pi \left(\begin{array}{c}\frac{1}{2}\\ -\frac{\sqrt{3}}{2}\end{array}\right).$$

$${\Omega }^{\ast }=\{{\theta }_1{\mathbf{k}}_1+{\theta }_2{\mathbf{k}}_2:-\frac{1}{2}\le {\theta }_j\le \frac{1}{2},j=1,2\},$$

$${\mathcal{L}}^A\Phi (\mathbf{x};\mathbf{k})$$

$$\Phi (\mathbf{x}+\mathbf{v};\mathbf{k})$$

$${\mathcal{L}}^A(\mathbf{k})\varphi (\mathbf{x};\mathbf{k})$$

$$\varphi (\mathbf{x}+\mathbf{v};\mathbf{k})$$

$$\mathcal{L}^{A}(\mathbf{k})=e^{-i\mathbf{k}\cdot\boldsymbol{x}}\mathcal{L}^{A}e^{i\mathbf{k}\cdot\boldsymbol{x}}=-(\nabla+i\mathbf{k})\cdot\big{(}A(\boldsymbol{x})(\nabla+i\mathbf{k})\big{)}.$$

$$0\le E_1(\mathbf{k})\le E_2(\mathbf{k})\le E_3(\mathbf{k})\le \cdots ,$$

$$E_b(\mathbf{k})\ge C>0,\text{for}|\mathbf{k}|\ge {\delta }_{b,1}\lambda .$$

$$\displaystyle f(\boldsymbol{x})=\frac{1}{|\Omega^{*}|}\sum_{b\geq 1}\int_{\Omega^{*}}\tilde{f}_{b}(\mathbf{k})\Phi_{b}(\boldsymbol{x};\mathbf{k})~{}d\mathbf{k},\quad\text{where}\quad\tilde{f}_{b}(\mathbf{k})=\big{\langle}\Phi_{b}(\cdot;\mathbf{k}),~{}f(\cdot)\big{\rangle},$$

$$\big{\|}f\big{\|}^{2}_{L^{2}(\mathbb{R}^{2})}=\frac{1}{|\Omega^{*}|}\sum_{b\geq 1}\int_{\Omega^{*}}|\tilde{f}_{b}(\mathbf{k})|^{2}~{}d\mathbf{k}~{}.$$

$$\Vert f{\Vert }_{H^s({\mathbb{R}}^2)}^2\approx$$

$$=$$

$$\approx$$

Full text

1 Introduction

This paper is concerned with a two dimensional (2D) wave equation

[TABLE]

where the material weight $W_{\varepsilon}$ has a slowly-modulated honeycomb structure,

[TABLE]

Here $A(\boldsymbol{x})$ defines a honeycomb structured media, $\varepsilon$ is a sufficiently small positive parameter, $\kappa(\varepsilon\boldsymbol{x})$ together with $B(\boldsymbol{x})$ are modulations applied to the honeycomb material, (see Section 2 for details). The 2D wave equation arises in describing the electromagnetic waves in an optical metamaterial whose permittivity and/or permeability are designed to have a hexagonal symmetry. These materials are also referred to as “photonic graphene”. In applications, modulations are often added to the honeycomb structure, for instance, domain-wall-like defects. These modulations leads to numerous topological phenomena [14, 24, 28, 30, 31, 32]

The subtle property of a honeycomb structured media defined by $A(\boldsymbol{x})$ is the existence of two-fold degenerate points, called Dirac points $(\mathbf{K}_{*},E_{D})$, lying in the dispersion surfaces of the elliptic operator \mathcal{L}^{A}\equiv-\nabla\cdot\big{(}A(\boldsymbol{x})\nabla\big{)}. The corresponding eigenfunctions, $\Phi_{1}(\boldsymbol{x})$ and $\Phi_{2}(\boldsymbol{x})=\overline{\Phi_{1}(-\boldsymbol{x})}$ satisfy

[TABLE]

Due to the conical property of Dirac points, the wave packets associated with such points have very novel propagation patterns. The goal of this work is to investigate the effective dynamics of such wave packets in (1.1) by deriving the envelope equation with a rigorous justification. To this end, we consider the wave equation (1.1) with the initial condition

[TABLE]

where $a_{j}(\cdot),~{}b_{j}(\cdot),~{}j=1,~{}2$ are all in Schwartz space $\mathcal{S}(\mathbb{R}^{2})$. This initial condition represents a general wave-packet spectrally localized at Dirac point $(\mathbf{K}_{*},E_{D})$. Thanks to the linearity of (1.1), the solution has two branches which evolve independently. Thus it is sufficient to consider the initial condition corresponding to one of the branches 111The analysis for the other branch is essentially the same by just changing $\sqrt{E_{D}}$ to $-\sqrt{E_{D}}$.

[TABLE]

where $\alpha_{10}(\cdot),~{}\alpha_{20}(\cdot)$ are in Schwartz space $\mathcal{S}(\mathbb{R}^{2})$. We remark that the factor $\varepsilon$ in front of $\psi(\boldsymbol{x},0)$ and $\partial_{t}\psi(\boldsymbol{x},0)$ are not essential due to the linearity of our problem. But this choice ensures that which brings great convenience to our analysis, i.e., $\|\psi(\boldsymbol{x},0)\|_{H^{s}(\mathbb{R}^{2})}=\mathcal{O}(1)$, $\|\partial_{t}\psi(\boldsymbol{x},0)\|_{H^{s}(\mathbb{R}^{2})}=\mathcal{O}(1)$.

We shall prove that the wave equation (1.1) with initial condition (1.5) has the following asymptotic solution,

[TABLE]

where $\alpha_{j}(\varepsilon\boldsymbol{x},\varepsilon t),~{}j=1,~{}2$ represent the slowly varying envelopes and $\eta(\boldsymbol{x},t)$ is the small correction to the leading order approximation. The main result of this paper stated in Theorem 4.1 is as follows. The envelope \alpha_{j}(\boldsymbol{X},T),~{}j=1,~{}2~{}\big{(}\boldsymbol{X}=(X_{1},~{}X_{2})=\varepsilon\boldsymbol{x},~{}T=\varepsilon t\big{)} satisfy the so-called Dirac equation with a varying mass

[TABLE]

with the initial condition

[TABLE]

Furthermore, given $\varepsilon>0$ sufficiently small, $\eta(\boldsymbol{x},t)$ can be controlled in Sobolev space $H^{s}(\mathbb{R}^{2})$ over a large but finite time scale $t\in[0,~{}\rho\varepsilon^{{}^{-1}}]$ as follows

[TABLE]

for any $0<\nu<1,~{}s\geq 0$ and $\rho>0$.

The main idea of our proof is inspired from the pioneer work by Fefferman and Weinstein [20]. They gave a rigorous justification of massless Dirac equation derived from the Schrödinger equation with a perfect honeycomb potential which corresponds to the special case $\kappa(\varepsilon\boldsymbol{x})\equiv 0$. In our current work, the case where $\kappa(\varepsilon\boldsymbol{x})\neq 0$ is considered in the wave equation (1.1) which has a second order derivative in time. These differences will bring technical difficulties in analysis which require new treatments.

In past decades of years, one of the most popular research subject is to understand and realize topological phenomena in different topological materials. One of the most successful example is the honeycomb-based material [22, 33, 38]. It stimulates the mathematical analysis of the Schrödinger equation with a honeycomb potential over the past few years [1, 2, 5, 9, 19, 20]. Fefferman and Weinstein rigorously proved the existence of Dirac points of Schrödinger operator with a generic honeycomb potential [19], and later gave a mathematical justification of the massless Dirac equation which governs the dynamics of the wave packets associated with Dirac points [20]. Topological edges states, strong binding limit, nonlinearity, and other aspects regarding this equation have also been investigated by them and others [3, 4, 5, 16, 17, 18, 21, 23, 27]. Meanwhile, many “artificial graphenes”, analogies of graphene in other fields, have been created to realize similar properties. Amongst those, photonic graphene has attracted a lot of interest due to its potential applications and relatively simple experimental realizations [24, 28, 32, 36]. To study electromagnetic waves in photonic graphene, we have to deal with Maxwell’s equations in [12, 13]. In a simple physical setting, for example, the propagation of transverse electronic fields can be described by the 2D wave equation (1.1). In a previous work, Lee-Thorp, Weinstein and Zhu [30] rigorously proved the existence of Dirac points of 2D elliptic operator with a honeycomb structured material weight and the existence of topological edge states in a domain-wall-modulated honeycomb structure. Their work paved the way to the mathematical analysis of electromagnetic waves in topological photonic materials. However, the wave packet dynamics has not been thoughtfully studied yet. In this present work, we shall give a rigorous investigation along this interesting and important line.

Before proceeding, we use the following basic requirements on the material weight.

Definition 1.1**.**

The $2\times 2$ complex-valued matrix function $W(\boldsymbol{x})$ is called a material weight, if

$W(\boldsymbol{x})$* is Hermitian and smooth for all $\boldsymbol{x}\in\mathbb{R}^{2}$,* 2. 2.

$W(\boldsymbol{x})$* is elliptic, i.e., for any $\boldsymbol{\xi}\in\mathbb{C}^{2}$, $\exists~{}0<C_{1}\leq C_{2}<+\infty$, such that $C_{1}|\boldsymbol{\xi}|^{2}\leq\bar{\boldsymbol{\xi}}\cdot W(\boldsymbol{x})\boldsymbol{\xi}\leq C_{2}|\boldsymbol{\xi}|^{2}$ for all $\boldsymbol{x}\in\mathbb{R}^{2}$.*

The rest of the paper is organized as follows. In section 2, we briefly review the basic Floquet-Bloch theory for 2D periodic elliptic operator, honeycomb structured media, Dirac points and other results which are used in later analysis; In section 3, we present the well-posedness of the envelope equation (1.7) in Schwartz space and its connection to topological edge states; In section 4, we conclude the main result$-$Theorem 4.1; In section 5, we give the detailed proofs of the key estimates which are essential to the proof of Theorem 4.1. In the appendix, we discuss how to apply our analysis to the non-modulating case, i.e., $\kappa(\varepsilon\boldsymbol{x})\equiv 0$, and give the detailed estimate of the solution to Dirac equation (1.7) in Schwartz space.

The following notations and conventions are used in this work:

$\Omega$ is the fundamental cell in (2.2), and $\Omega^{*}$ the dual fundamental cell defined in (2.6). 2. 2.

The operators $\mathcal{L}^{A},~{}\mathcal{L}^{B},~{}\mathcal{L}^{B}_{\kappa}$ are denoted as follows:

[TABLE] 3. 3.

The standard notations for function spaces are used:

• (1)

$H^{s}(\mathbb{R}^{2})$ is the Sobolev space, i.e., $f(\boldsymbol{x})\in H^{s}(\mathbb{R}^{2})$, $\|f\|_{H^{s}(\mathbb{R}^{2})}^{2}=\sum\limits_{|\mathbf{n}|\leq s}\int_{\mathbb{R}^{2}}|\partial^{\mathbf{n}}_{\boldsymbol{x}}f(\boldsymbol{x})|^{2}d\boldsymbol{x}<\infty$;

• (2)

$\mathcal{S}(\mathbb{R}^{2})$ is the Schwartz space, i.e., $f(\boldsymbol{x})\in\mathcal{S}(\mathbb{R}^{2})$, $\sup\limits_{\boldsymbol{x}\in\mathbb{R}^{2}}|\boldsymbol{x}^{\mathbf{m}}\partial_{\boldsymbol{x}}^{\mathbf{n}}f(\boldsymbol{x})|<\infty,~{}\forall~{}\mathbf{m},\mathbf{n}\in\mathbb{N}^{2}$;

• (3)

$C_{b}^{\infty}(\mathbb{R}^{2})$ contains all functions which are bounded and smooth. 4. 4.

We use the indicator function $\chi(D)$ to define the Bloch spectral cutoff, i.e., for any constant $C\geq 0$, and $f:\Omega^{*}\to\mathbb{C}$,

[TABLE] 5. 5.

The Pauli matrices $\sigma_{j},~{}j=1,~{}2,~{}3$ are:

[TABLE] 6. 6.

$A\approx B$ if and only if there exist two constants $0<C_{1}\leq C_{2}$ such that $C_{1}A\leq B\leq C_{2}A$. 7. 7.

We use notations \big{\langle}\cdot,\cdot\big{\rangle} and \big{\langle}\cdot,\cdot\big{\rangle}_{L^{2}(\Omega)} to distinguish the inner products on $L^{2}(\mathbb{R}^{2})$ and $L^{2}(\Omega)$, i.e.,

[TABLE] 8. 8.

For any $f(\boldsymbol{x})\in L^{2}(\mathbb{R}^{2})$, $\hat{f}(\boldsymbol{\xi})$ represents its Fourier transform while $\tilde{f}_{b}(\mathbf{k})$ stands for its Bloch component, i.e.,

[TABLE] 9. 9.

The repeated index summation convention is used throughout.

2 Preliminaries

In this section, we list the Floquet-Bloch theory, honeycomb structure, Dirac points which are desired for the arguments of this work. We refer readers to [8, 10, 11, 15, 29, 30, 34, 39] for more details.

2.1 Floquet-Bloch theory

A lattice $\Lambda$ in $\mathbb{R}^{2}$ is generated by two linearly independent vectors $\mathbf{v}_{1}$, $\mathbf{v}_{2}$, i.e.,

[TABLE]

and the fundamental cell is chosen to be the parallelogram:

[TABLE]

where $|\Omega|$ is the area of $\Omega$. In this work, we specify the lattice $\Lambda$ to be a triangular lattice with lattice vectors

[TABLE]

The dual lattice

[TABLE]

is generated by the dual lattice vectors $\mathbf{k}_{1},~{}\mathbf{k}_{2}$ which satisfy $\mathbf{k}_{i}\cdot\mathbf{v}_{j}=2\pi\delta_{ij},~{}i,~{}j=1,~{}2$. For the triangular lattice defined by (2.3), the dual lattice vectors are

[TABLE]

Throughout this work, we choose the parallelogram $\Omega^{*}$:

[TABLE]

as the fundamental dual cell.222In many literatures, another choice of the fundamental cell is the Brillouin zone $\mathcal{B}$, consisting the points $\mathbf{k}\in\mathbb{R}^{2}$ which are closer to the origin than to any other lattice points in $\Lambda^{*}$. For the triangular lattice, the Brillouin zone is a hexagon, see for example [19].

Let $L^{2}_{\text{per}}(\mathbb{R}^{2}/\Lambda)$ denote a subspace of $L^{2}_{\text{loc}}(\mathbb{R}^{2})$ containing all $\Lambda$-periodic functions, namely $f\in L^{2}_{\text{loc}}(\mathbb{R}^{2})$ and $f(\boldsymbol{x}+\mathbf{v})=f(\boldsymbol{x})$, $\forall~{}\boldsymbol{x}\in\mathbb{R}^{2},~{}\forall~{}\mathbf{v}\in\Lambda$. For each $\mathbf{k}\in\mathbb{R}^{2}$, we denote $f\in L_{\mathbf{k}}^{2}(\mathbb{R}^{2}/\Lambda)$ if $e^{-i\mathbf{k}\cdot\boldsymbol{x}}f(\boldsymbol{x})\in L^{2}_{\text{per}}(\mathbb{R}^{2}/\Lambda),~{}i.e.,~{}f(\boldsymbol{x}+\mathbf{v})=e^{i\mathbf{k}\cdot\mathbf{v}}f(\boldsymbol{x})$, for all $\boldsymbol{x}\in\mathbb{R}^{2},~{}\mathbf{v}\in\Lambda$. Similarly, we can define the Sobolev space $H^{s}_{\mathbf{k}}(\mathbb{R}^{2}/\Lambda),~{}s\geq 0$.

Suppose the matrix function $A(\boldsymbol{x})$ is a material weight in the sense of Definition 1.1, and further $\Lambda$-periodic. \mathcal{L}^{A}=-\nabla\cdot\big{(}A(\boldsymbol{x})\nabla\big{)} is an elliptic operator with periodic coefficient, thus the Floquet-Bloch theory applies. For any $\mathbf{k}\in\mathbb{R}^{2}$, consider the following $L^{2}_{\mathbf{k}}$-Floquet-Bloch elliptic eigenvalue problem with pseudo-periodic boundary condition

[TABLE]

Since the above eigenvalue problem is invariant under the translation $\mathbf{k}\rightarrow\mathbf{k}+\mathbf{k}^{\prime}$, $\forall~{}\mathbf{k}^{\prime}\in\Lambda^{*}$, it is sufficient to just pay attention to $\mathbf{k}$ varying over $\Omega^{*}$. Alternatively, one can obtain the periodic elliptic boundary problem by setting $\Phi(\boldsymbol{x};\mathbf{k})=e^{i\mathbf{k}\cdot\boldsymbol{x}}\phi(\boldsymbol{x};\mathbf{k}),~{}\mathbf{k}\in\Omega^{*}$, i.e.,

[TABLE]

where

[TABLE]

Notice that $\mathcal{L}^{A}(\mathbf{k})$ is a self-adjoint elliptic operator in $L^{2}_{\text{per}}(\mathbb{R}^{2}/\Lambda)$. For each fixed $\mathbf{k}\in\Omega^{*}$, the above eigenvalue problem has a series of discrete spectrum (or eigenvalues) [8, 20, 39]:

[TABLE]

with eigenpairs \big{(}\phi_{b}(\boldsymbol{x};\mathbf{k}),~{}E_{b}(\mathbf{k})\big{)},~{}b\geq 1, where $\{\phi_{b}(\boldsymbol{x};\mathbf{k})\}_{b\geq 1}$ can be chosen a complete orthogonal basis in $L^{2}_{\text{per}}(\mathbb{R}^{2}/\Lambda)$. The eigenvalues $E_{b}(\mathbf{k})$ are called band dispersion functions or Bloch bands which are Lipschitz continuous, and $E_{b}(\mathbf{k})=0$ if and only if $b=1$ and $\mathbf{k}=\boldsymbol{0}$ with the corresponding normalized eigenfunction $\Phi_{1}(\boldsymbol{x};\boldsymbol{0})=\phi_{1}(\boldsymbol{x};\boldsymbol{0})\equiv|\Omega|^{-\frac{1}{2}}$. Thus, for any $\lambda>0$, there exist a positive constant $C$ such that

[TABLE]

We will choose an appropriate constant $\lambda$ for the proof convenience in Section 5.

The corresponding quasi-periodic eigenfunctions $\Phi_{b}(\boldsymbol{x};\mathbf{k})$ are called Bloch modes. For any given $\mathbf{k}\in\Omega^{*}$, $\Phi_{b}(\boldsymbol{x};\mathbf{k})~{}(\text{or}~{}\phi_{b}(\boldsymbol{x};\mathbf{k})),~{}b\geq 1$ is analytic for $\boldsymbol{x}\in\mathbb{R}^{2}$ from the regularity theory of elliptic operator. Moreover, the set of all Bloch modes $\{\Phi_{b}(\boldsymbol{x};\mathbf{k})\}_{b\geq 1,~{}\mathbf{k}\in\Omega^{*}}$ is complete in $L^{2}(\mathbb{R}^{2})$. That is for any $f\in L^{2}(\mathbb{R}^{2})$,

[TABLE]

and the Parseval formula for the Bloch decomposition in $L^{2}(\mathbb{R}^{2})$ holds,

[TABLE]

Thanks to the Weyl’s law, i.e., $E_{b}(\mathbf{k})\approx b~{}(b\gg 1)$ uniformly for all $\mathbf{k}\in\Omega^{*}$, and then any $f\in H^{s}(\mathbb{R}^{2})$, $\|f\|_{H^{s}(\mathbb{R}^{2})}$ can be approximated by

[TABLE]

2.2 Honeycomb structured material weight and Dirac points

Before introducing the honeycomb structured material weight, we define the following symmetry operators acting on a function $f(\boldsymbol{x})$ defined in $\mathbb{R}^{2}$. Parity inversion operator $\mathcal{P}$: $\mathcal{P}[f](\boldsymbol{x})=f(-\boldsymbol{x})$; Complex conjugate operator $\mathcal{C}$: $\mathcal{C}[f](\boldsymbol{x})=\overline{f(\boldsymbol{x})}$; $\frac{2\pi}{3}-$rotation operator $\mathcal{R}$: $\mathcal{R}[f](\boldsymbol{x})=f(R^{*}\boldsymbol{x})$ where $R$ is $\frac{2\pi}{3}-$clockwise rotation matrix

[TABLE]

A honeycomb structured material weight is defined as follows

Definition 2.1**.**

A $2\times 2$ complex-valued matrix function $A(\boldsymbol{x})$ is called a honeycomb structured material weight if it satisfies the Definition 1.1 and further

$A(\boldsymbol{x})$* is $\Lambda$-periodic, i.e., $A(\boldsymbol{x}+\mathbf{v})=A(\boldsymbol{x}),~{}\forall~{}\boldsymbol{x}\in\mathbb{R}^{2},~{}\mathbf{k}\in\Lambda$;* 2. 2.

$A(\boldsymbol{x})$* is $\mathcal{PC}$-invariant, i.e., $\mathcal{PC}[A](\boldsymbol{x})\equiv\overline{A(-\boldsymbol{x})}=A(\boldsymbol{x}),~{}\forall~{}\boldsymbol{x}\in\mathbb{R}^{2}$;* 3. 3.

$A(\boldsymbol{x})$* satisfies $\mathcal{R}[A](\boldsymbol{x})\equiv A(R^{*}\boldsymbol{x})=R^{*}A(\boldsymbol{x})R,~{}\forall~{}\boldsymbol{x}\in\mathbb{R}^{2}$.*

A consequence of $A(\boldsymbol{x})$ being a honeycomb structured material weight is that $\mathcal{PC}\mathcal{L}^{A}=\mathcal{L}^{A}\mathcal{PC},~{}\mathcal{R}\mathcal{L}^{A}=\mathcal{L}^{A}\mathcal{R}$. Generically, this leads to the existence of Dirac points in the dispersion surfaces of \mathcal{L}^{A}=-\nabla\cdot\big{(}A(\boldsymbol{x})\nabla\big{)}. The specific definition of Dirac points is given as follows, see [19, 30].

Definition 2.2**.**

The quasi-momentum/eigenvalue pair $(\mathbf{K}_{*},E_{D})\in\Omega^{*}\times\mathbb{R}_{+}$ is called a Dirac point if there exists an integer $b_{*}\geq 1$ and Floquet-Bloch eigenpairs mappings

[TABLE]

such that:

1. (1)

$E_{D}=E_{b_{*}}(\mathbf{K}_{*})=E_{b_{*}+1}(\mathbf{K}_{*})$* is a two-fold degenerate $L^{2}_{\mathbf{K}_{*}}$– eigenvalue of $\mathcal{L}^{A}$. There exists two orthogonal eigenfunctions $\Phi_{1}(\boldsymbol{x}),~{}\Phi_{2}(\boldsymbol{x})=\overline{\Phi_{1}(-\boldsymbol{x})}$,*

[TABLE] 2. (2)

Denote that $E_{-}(\mathbf{k})=E_{b_{*}}(\mathbf{k}),~{}E_{+}(\mathbf{k})=E_{b_{*}+1}(\mathbf{k})$, and $\boldsymbol{\kappa}=(\kappa_{1},\kappa_{2})=\mathbf{k}-\mathbf{K}_{*}$. There exist $v_{{}_{F}}>0$ and $q_{0}>0$, such that for $0\leq|\boldsymbol{\kappa}|\leq q_{0}$,

[TABLE]

where $|e_{\pm}(\boldsymbol{\kappa})|\leq C|\boldsymbol{\kappa}|$ for some constant $C$.

Consider the two high symmetric points in $\Omega^{*}$, $\mathbf{K}=\frac{1}{3}(\mathbf{k}_{1}-\mathbf{k}_{2})$ and $\mathbf{K}^{\prime}=-\mathbf{K}$. In the previous literature, Lee-Thorp, Weinstein and Zhu proved, if $A(\boldsymbol{x})$ is a honeycomb structured material weight in the sense of Definition 2.1, the Dirac points generically appear in two adjacent dispersion surfaces of $\mathcal{L}^{A}$ conically intersecting at $\mathbf{K}_{*}=\mathbf{K}$ and $\mathbf{K}^{\prime}$, see Theorem 4.2 and 4.10 in [30].

Furthermore, we have two more conclusions at the following proposition

Proposition 2.3**.**

Let $(\mathbf{K}_{*},E_{D})$ be a Dirac point in the sense of Definition 2.2. If $b\notin\{+,-\}$, there exists a constant $C>0$, and $q_{1}>0$ small, such that for any $0\leq|\mathbf{k}-\mathbf{K}_{*}|\leq q_{1}$

[TABLE]

Let $\Phi_{-}(\boldsymbol{x};\mathbf{k})=\Phi_{b_{*}}(\boldsymbol{x};\mathbf{k}),~{}\Phi_{+}(\boldsymbol{x};\mathbf{k})=\Phi_{b_{*}+1}(\boldsymbol{x};\mathbf{k})$. When $0<|\boldsymbol{\kappa}|\leq q_{0}$, one can expand $\Phi_{\pm}(\boldsymbol{x};\mathbf{k})$ in the following form,

[TABLE]

The proof of (2.19) in Proposition 2.3 is a direct consequence of (2.11) and the Lipschitz continuity for all eigenvalues, while for the rigorous proof of (2.20) is referred to Theorem 3.2 in [20].

In this work, we consider a slowly modulated honeycomb structured material weight $W_{\varepsilon}(\boldsymbol{x})=A(\boldsymbol{x})+\varepsilon\kappa(\varepsilon\boldsymbol{x})B(\boldsymbol{x})$. Throughout this work, we require the following assumptions on $W_{\varepsilon}(\boldsymbol{x})$ hold

Assumption 1**.**

The $2\times 2$ matrix $W_{\varepsilon}(\boldsymbol{x})=A(\boldsymbol{x})+\varepsilon\kappa(\varepsilon\boldsymbol{x})B(\boldsymbol{x})$ is a material weight in the sense of Definition 1.1, and

1. ($A_{1}$)

$A(\boldsymbol{x})$* is a honeycomb structured material weight in the sense of Definition 2.1;* 2. ($A_{2}$)

$\kappa(\boldsymbol{X})$* is a real scalar function in $C^{\infty}_{b}(\mathbb{R}^{2})$;* 3. ($A_{3}$)

$B(\boldsymbol{x})$* is smooth, Hermitian, $\Lambda$-periodic, and $B(-\boldsymbol{x})=B(\boldsymbol{x})$, $\overline{B(\boldsymbol{x})}=-B(\boldsymbol{x})$;*

Remark: Assumption $(A_{3})$ implies that $\mathcal{PC}\mathcal{L}^{B}=-\mathcal{L}^{B}\mathcal{PC}$ and further $B(\boldsymbol{x})=b(\boldsymbol{x})\sigma_{2}$ with $b(\boldsymbol{x})$ real and even. This assumption is related to the time-reversal symmetry breaking in the study of photonic topological insulators in real applications [24, 28, 30]. It also indicates that the second order operator $\mathcal{L}^{B}_{\kappa}$ is actually a first order operator on $\eta$, i.e.,

[TABLE]

In practical applications, there is another way to break $\mathcal{PC}$-symmetry by assuming $B(\boldsymbol{x})=\tilde{b}(\boldsymbol{x})I_{2\times 2}$ with $\tilde{b}(\boldsymbol{x})$ real and odd, which is related to parity symmetry breaking [30]. In this case $\mathcal{L}^{B}_{\kappa}$ remains a second order operator and our analysis does not apply.

A typical example of a media satisfying Assumption 1 in real applications is the magneto-optical material given in Haldane and Ragu’s work [24]. In their physical setup, the material weights written in our notation are

[TABLE]

where $\epsilon(\boldsymbol{x})$ is the electric permittivity which is even, real and $\mathcal{R}$-invariant, and $\gamma(\varepsilon\boldsymbol{x})$ represents the Fardy-rotation effect which is assumed to be small and slowly varying.

We end this section by including the following important results [30].

Proposition 2.4**.**

Suppose that $\Phi_{1}(\boldsymbol{x}),~{}\Phi_{2}(\boldsymbol{x})=\overline{\Phi_{1}(-\boldsymbol{x})}$ are eigenfunctions of $\mathcal{L}^{A}$ with respect to the Dirac point $(E_{D},\mathbf{K}_{*})$ given in Definition 2.2, $B(\boldsymbol{x})$ satisfies $(A_{3})$, and define

[TABLE]

Then, the following identities hold:

[TABLE]

and

[TABLE]

Here we define \vartheta_{\sharp}=\big{\langle}\Phi_{1},~{}\mathcal{L}^{B}\Phi_{1}\big{\rangle}_{L^{2}(\Omega)} and assume $\vartheta_{\sharp}\neq 0$ in this work.

3 Dirac equation with a varying mass

This paper is to show that the Dirac equation (1.7) governs the envelope dynamics of the wave problem (1.1) under a prescribed initial condition (1.5). In this section, we first state results on well-posedness, estimates on solutions to the Dirac equation (1.7) in the following proposition.

Proposition 3.1**.**

Let $E_{D}>0$, $v_{{}_{F}}>0$, $\vartheta_{\sharp}\neq 0$ be given constants as before, $\kappa(\boldsymbol{X})\in C_{b}^{\infty}(\mathbb{R}^{2})$, and \boldsymbol{\alpha}_{0}(\boldsymbol{X})=\big{(}\alpha_{10}(\boldsymbol{X}),\alpha_{20}(\boldsymbol{X})\big{)}\in\mathcal{S}(\mathbb{R}^{2}). Then, for any $0<\rho<+\infty,~{}s>2$, the Dirac equation (1.7) has a unique solution \boldsymbol{\alpha}(\boldsymbol{X},T)=\big{(}\alpha_{1}(\boldsymbol{X},T),~{}\alpha_{2}(\boldsymbol{X},T)\big{)}^{T}\in C^{\infty}\big{(}[0,\infty)\times\mathbb{R}^{2}\big{)} and

[TABLE]

Moreover, for any $T\geq 0,~{}l\in\mathbb{N}$, $\partial_{T}^{l}\alpha_{j}(\cdot,T)\in\mathcal{S}(\mathbb{R}^{2})$, specifically, $\forall~{}M\in\mathbb{N}$, $\mathbf{n}\in\mathbb{N}^{2}$, there exists a constant $C>0$ such that

[TABLE]

Note that Dirac equation (1.7) is actually a first order linear hyperbolic system. If $\kappa(\boldsymbol{X})$ is a constant, (1.7) has a unique solution in Schwartz space for the time $T\in[0,+\infty)$ in term of the Fourier transform arguments [20]. However, for a general $\kappa(\boldsymbol{X})\in C_{b}^{\infty}(\mathbb{R}^{2})$, a comprehensive proof to the above proposition will be postponed in the Appendix B.

One of the most interesting applications of the reduced Dirac equation (1.7) is its capability to describe dynamics of topological edge states. To be more specific, suppose that $\kappa(\boldsymbol{X})$ is a domain wall function, i.e., $\kappa(\boldsymbol{X})=\tilde{\kappa}(\zeta)$ with $\tilde{\kappa}(\zeta)\to\pm\kappa_{\infty}$ as $\zeta\to\pm\infty$, where $\zeta=\boldsymbol{\mathfrak{K}}\cdot\boldsymbol{X}\in\mathbb{R}$ and $\boldsymbol{\mathfrak{K}}=(\mathfrak{K}_{1},\mathfrak{K}_{2})\in\Lambda^{*}$, $\kappa_{\infty}>0$, see [17, 30] for details. Hereafter, we drop the tilde on top of $\kappa$.

Let

[TABLE]

Rewriting the Dirac equation (1.7) in new coordinates system $(\xi,\zeta)$ yields that

[TABLE]

where $\alpha_{j}=\tilde{\alpha}_{j}(\xi,\zeta,T)=\alpha_{j}(\boldsymbol{x}(\xi,\zeta),T)$ and we have used the same symbols before and after changing coordinates for notational convenience.

We are interested in a particular solution to (3.3) which decays to zero as $|\zeta|\to\infty$ and keeps periodic in $\xi-$direction. This solution is referred to as the topological edge state. Namely, let

[TABLE]

It deduces that $\boldsymbol{\beta}(\zeta;k_{\parallel})$ satisfies the following eigenvalue problem in $L^{2}(\mathbb{R})$,

[TABLE]

where $\mathcal{D}(k_{\parallel})$ is the 1D Dirac operator

[TABLE]

In the previous work, Lee-Thorp, Weinstein and Zhu [30] derived the same equation as (3.4) in the case that $k_{\parallel}=0$. In this scenario, there exist the so-called zero-energy states for $\mathcal{D}(0)$. That is, \big{(}\boldsymbol{\beta}(\zeta;0),0\big{)} is the eigenpair of the operator $\mathcal{D}(0)$ in $L^{2}(\mathbb{R})$ with

[TABLE]

and here $\gamma$ is the normalization constant.

Our derivation and analysis demonstrate the existence of edge states proved in [30] from the evolutionary effective envelope equation (1.7). Equation (1.7) can be further used to study the dynamics of such states as well as their interactions with defects, perturbations by manipulating $\kappa(\zeta,\xi)$.

In order to show that the envelope equation (1.7) can exhibit many interesting solutions which describe the novel and subtle physical phenomena, we illustrate a typical propagation pattern by solving (1.7) numerically. In this simulation, the coefficients of the equation is normalized for simplicity. Alternatively, we set $\frac{v_{{}_{F}}}{2\sqrt{E_{D}}}=1,~{}\frac{\vartheta_{\sharp}}{2\sqrt{E_{D}}}=1$. The modulation that we choose is

[TABLE]

The initial condition is

[TABLE]

In Figure 1, we plot several snapshots of the intensity of the solution to equation (1.7) at three successive times. It can be seen that the waves travel along the edge without any energy leaking to the bulk or traveling back. This interesting phenomenon is related to the topologically protected wave propagation which is one of the current focuses in many applied fields. Our rigorous justification of (1.7) from (1.1) provides a solid mathematical foundation for such interesting problems. Due to the length and scope of this paper, we leave the further analysis on the reduced envelope equation (1.7) in future works.

4 Main results

The main goal of this paper is to show that the equation (1.1) with initial condition (1.5) has the asymptotic solution (1.6) with the envelopes $\alpha_{j}(\varepsilon\boldsymbol{x},\varepsilon t),~{}j=1,~{}2$ satisfying the Dirac equation (1.7)(1.8). The well-posedness of Cauchy problem (1.1)(1.5) is a standard result by the theory on linear hyperbolic systems, see e.g., [26, 35]. Our task is reduced to rigorously justify that the error $\eta(\boldsymbol{x},t)$ is small over a large but finite time. To this end, we substitute (1.6) into (1.1)(1.5) and obtain the equation of $\eta(\boldsymbol{x},t)$,

[TABLE]

with the initial condition

[TABLE]

where $F_{1}(\boldsymbol{x},t)$ and $F_{2}(\boldsymbol{x},t)$ are

[TABLE]

[TABLE]

The main result is concluded as follows.

Theorem 4.1**.**

Suppose $W_{\varepsilon}(\boldsymbol{x})$ satisfies Assumption 1, the operator \mathcal{L}^{A}=-\nabla\cdot\big{(}A(\boldsymbol{x})\nabla\big{)} has a Dirac point $(\mathbf{K}_{*},E_{D})$, $\Phi_{1}(\boldsymbol{x})$ and $\Phi_{2}(\boldsymbol{x})$ are the associated eigenfunctions given in Definition 2.2, and $\alpha_{10}(\boldsymbol{X})$, $\alpha_{20}(\boldsymbol{X})$ are Schwartz functions of $\boldsymbol{X}\in\mathbb{R}^{2}$. Then the wave equation (1.1) with the initial condition (1.5) has a unique solution of the form (1.6), where $\alpha_{j}(\boldsymbol{X},T),~{}j=1,~{}2$ are the solution of the system (1.7)(1.8), and for any $s\geq 0$, $\rho>0$, $0<\nu<1$ and $\varepsilon>0$ sufficiently small,

[TABLE]

Here $C$ is independent of $\varepsilon$.

Remark: In the literatures on wave packet problems, the macroscopic reference frame $\boldsymbol{X}=\varepsilon\boldsymbol{x},~{}T=\varepsilon t$ with a $\varepsilon-$scaled Sobolev space $H^{s}_{\varepsilon}(\mathbb{R}^{2})$, i.e., $\forall~{}f(\boldsymbol{X})\in H^{s}_{\varepsilon}(\mathbb{R}^{2})$,

[TABLE]

is frequently used, see for instance [5, 23, 37]. As shown in Theorem 4.1, we do the error estimates in the microscopic scales with the regular Sobolev spaces. Actually, the two treatments are essentially equivalent. Indeed, due to the linearity, the equation remain the same under rescaling $\widetilde{\psi}=\frac{\psi}{\varepsilon}$. The following identity shows the equivalence

[TABLE]

Proof of Theorem 4.1. We follow the standard procedure for wave packet problems. Namely, we first spectrally decompose the error by the Floquet-Bloch theory. Then the spectral components are estimated separately.

Recalling the completeness of Bloch modes of $\mathcal{L}^{A}$ in $L^{2}(\mathbb{R}^{2})$, we have

[TABLE]

where the error component

[TABLE]

Then, for any $b\geq 1$, $\tilde{\eta}_{b}(\mathbf{k},t)$ satisfies

[TABLE]

with initial condition

[TABLE]

By Duhamel’s principle, we rewrite (4.8) as the integral form

[TABLE]

Especially, when $b=1$ and $\mathbf{k}=\boldsymbol{0}$,

[TABLE]

the fourth term on the right hand side vanishes since $\Phi_{1}(\boldsymbol{x};\boldsymbol{0})$ is a constant.

From (4)-(4), we can see that the error components $\tilde{\eta}_{b}(\mathbf{k},t)$ appear in a very different way from those in [20]. First, $\frac{1}{\sqrt{E_{b}(\mathbf{k})}}$ appears in the error equation which brings secular terms at $b=1$ and $\mathbf{k}=\boldsymbol{0}$ as shown in (4). The other main difference is the presence of implicit terms $\tilde{\mathcal{Q}}_{b}[\eta](\mathbf{k},t)$ which are caused by the effects of modulation/perturbation. Consequently, more efforts and techniques are desired. In our analysis below, we deal with the implicit terms by Gronwall’s inequality after proving the boundedness of the operator $\frac{\mathcal{L}^{B}_{\kappa}}{\sqrt{\mathcal{L}^{A}}}$ in $H^{s}(\mathbb{R}^{2})$. To handle the singularity, we carefully estimate the error components near and away from the singularity separately.

Note that if $\kappa(\varepsilon\boldsymbol{x})\equiv 0$, which happened similarly in Fefferman and Weinstein’s work [20] on the Schrödinger equation, then $\tilde{\mathcal{Q}}_{b}[\eta](\mathbf{k},t),~{}b\geq 1$ vanished. As a consequence, each $\tilde{\eta}_{b}(\mathbf{k},t)$ is explicitly represented without coupling to other components.

In general, (4) and (4) imply that $\eta(\boldsymbol{x},t)$ satisfies the integral equation

[TABLE]

and more precisely, $G_{j}(\boldsymbol{x},t),~{}j=0,~{}1,~{}2$ and $\mathcal{Q}[\eta](\boldsymbol{x},t)$ satisfy

[TABLE]

To achieve the error bound of $\eta(\boldsymbol{x},t)$ from (4), we require the following two important propositions.

Proposition 4.2**.**

According to (4)-(4.13) and the Floquet-Bloch theory (2.13), for any $t\geq 0$ and given the integer $s\geq 0$, one can get

[TABLE]

Proof. We first recall the well-known results on the Riesz transform \nabla\big{(}\mathcal{L}^{A}\big{)}^{-\frac{1}{2}} which ensures the boundedness of the operator \nabla\big{(}\mathcal{L}^{A}\big{)}^{-\frac{1}{2}} from $L^{2}(\mathbb{R}^{2})$ to $L^{2}(\mathbb{R}^{2})$, see [6, 7, 25] and the references therein for details. Since the operator $\mathcal{L}^{A}$ is self-adjoint, for any two functions $f\in H^{1}(\mathbb{R}^{2})$, $\varphi\in\mathcal{S}(\mathbb{R}^{2})$ with $\|\varphi\|_{L^{2}(\mathbb{R}^{2})}=1$, one can obtain the following estimate by Riesz transform in the dual case,

[TABLE]

Then, for any $t\geq 0$, we take the absolute value of $\tilde{\mathcal{Q}}_{b}[\eta](\mathbf{k},t)$,

[TABLE]

By the Floquet-Bloch theory and Minkowski’s integral inequality, it follows that

[TABLE]

Further, we obtain the estimate in $H^{s}(\mathbb{R}^{2})$ for any integer $s\geq 1$,

[TABLE]

Now we turn to estimates of the first three terms in (4) which we conclude in the following Proposition.

Proposition 4.3**.**

As $G_{j}(\boldsymbol{x},t),~{}j=0,~{}1,~{}2$ are defined in (4.13), for any $s\geq 0$, $\rho>0$, $0<\nu<1$ and $\varepsilon>0$ sufficiently small, the following three statements hold

[TABLE]

where each constant $C$ is independent of $\varepsilon$.

We shall prove this proposition in the subsequential section.

Therefore, thanks to the above Proposition 4.2-4.3, we can now prove Theorem 4.1. For any $s\geq 0$, $0\leq t\leq\rho\varepsilon^{{}^{-1}}$ and $0<\nu<1$,

[TABLE]

By utilizing Gronwall’s inequality, we obtain

[TABLE]

and finally

[TABLE]

To complete the proof of Theorem 4.1, we only need to justify (4.16)-(4.18) of which the detailed verifications will be given successively in the next section.

5 Proof of Proposition 4.3

In this section, we shall give the detailed proof of Proposition 4.3. Hereafter, we just suppress the subscript of $\mathbf{K}_{*}$ as $\mathbf{K}$ for simplicity. For the convenience of our proof, we set $\lambda=\frac{|\mathbf{K}|}{2}$ which ensures the lower positive bound in (2.12), and let $0<q_{1}<|\mathbf{K}|-\lambda$ in the Proposition 2.3. Before proceeding further, we present several results concluded in Proposition 5.1 which will be frequently used in our proof.

Proposition 5.1**.**

Let $\Gamma(\boldsymbol{X})\in\mathcal{S}(\mathbb{R}^{2})$ and $\Psi(\boldsymbol{x})\in C^{\infty}(\mathbb{R}^{2})\cap L_{\mathbf{K}}^{2}(\mathbb{R}^{2}/\Lambda)$. Then, the following statements hold

[TABLE]

and

[TABLE]

Further, for any $\mathbf{k}\in\Omega^{*}$, $\frac{1}{2}<\nu_{1}<1$ and $\varepsilon>0$ sufficiently small,

[TABLE]

where each $C$ is a generic constant.

The detailed proof is omitted, we refer the readers to [20] for complete discussions.

Utilizing the Poisson-Summation (5.1) and the property (5.2) in Proposition 5.1, we can obtain the following lemma.

Lemma 5.2**.**

Let $\tilde{G_{j}}_{1}(\mathbf{k},t),~{}j=0,~{}1,~{}2$ be defined in (4). For any $\mathbf{k}\in\Omega^{*}$, $|\mathbf{k}|<\frac{|\mathbf{K}|}{2}$, $N>2$ and $0\leq t\leq\rho\varepsilon^{{}^{-1}}$,

[TABLE]

Proof. Without loss of generality, we only consider the case of $j=1$, the other two cases $j=0,~{}2$ can be treated similarly. For the convenience, we have to introduce the following notations,

[TABLE]

Thanks to Proposition 3.1, we can conclude that $\Gamma_{1r}=\Gamma_{1r}(\boldsymbol{X},T)\in\mathcal{S}(\mathbb{R}^{2})$ for any $T\geq 0$, and $\Psi_{1r}=\Psi_{1r}(\boldsymbol{x})\in C^{\infty}(\mathbb{R}^{2})\cap L^{2}_{\mathbf{K}}(\mathbb{R}^{2}/\Lambda),~{}r\in\{1,\cdots,6\}$. Then, Poisson-Summation (5.1) yields that

[TABLE]

If $|\mathbf{k}|<\frac{|\mathbf{K}|}{2}$, i.e., $|\mathbf{k}-\mathbf{K}|\geq q_{1}$, there exists a constant $C>0$ by (5.3) such that

[TABLE]

According to (5.2),

[TABLE]

Since $F_{1}(\boldsymbol{x},\tau)=\varepsilon^{2}\sum\limits_{r=1}^{6}\Gamma_{1r}(\varepsilon\boldsymbol{x},\varepsilon\tau)\cdot\Psi_{1r}(\boldsymbol{x})$, we get when $N>2$, $t\in[0,\rho\varepsilon^{{}^{-1}}]$,

[TABLE]

In the following justifications, we will frequently use the consequence in Proposition 3.1 and 5.1.

5.1 Proof of (4.16)

Let $N>2$, then by Lemma 5.2,

[TABLE]

(2.12), (4)-(4.13) imply

[TABLE]

Therefore, for any $t\geq 0$, we acquire that

[TABLE]

Thus (4.16) is justified.

5.2 Proof of (4.18)

As in the case of the proof above, it is straightforward to show that

[TABLE]

Moreover, in the calculations below we will use the fact that for any $0\leq\tau\leq\rho\varepsilon^{{}^{-1}}$, $s\geq 0$,

[TABLE]

Then, thanks to Lemma 5.2, we get for $0\leq t\leq\rho\varepsilon^{{}^{-1}}$,

[TABLE]

This leads to the result (4.18), i.e.,

[TABLE]

5.3 Proof of (4.17)

In this subsection, we turn to the key estimate (4.17). The main idea of spectral domain decomposition is similar to that presented in [20] except several considerable modifications due to the differences of the underlying problems. We first divide the Bloch components of $G_{1}(\boldsymbol{x},t)$ as follows,

[TABLE]

Since we have assumed that $\varepsilon>0$ sufficiently small satisfying $\varepsilon^{\frac{1}{2}}<\frac{|\mathbf{K}|}{2}$, and $0<q_{1}<\frac{|\mathbf{K}|}{2}$ as mentioned before, the second part on the right hand side of above can be divided into $b\in\{+,-\}$ part ${G_{1}}_{\mathrm{D}}$ and $b\notin\{+,-\}$ part ${G_{1}}_{\mathrm{D^{C}}}$. Further, we decompose ${G_{1}}_{\mathrm{D}}$ and ${G_{1}}_{\mathrm{D^{C}}}$ into their quasi-momentum components near and far away from $\mathbf{K}$. Specifically, we have the following decomposition,

[TABLE]

and

[TABLE]

Here $\chi$ is the indicator function, $\frac{1}{2}<\nu_{1}<1$ is to be determined and $q_{1}$ is a specified constant independent of $\varepsilon$ such that (2.19) holds.

Then by definition (2.14) and (2.15), $G_{1}(\boldsymbol{x},t)$ can be approximated as follows:

[TABLE]

The estimate of the first part on the right hand side of above directly follows from Lemma 5.2 when $N\geq 3$, and next we verify the error bounds of the other four parts term by term.

Recall that $\tilde{G_{1}}_{b}(\mathbf{k},t)$ is defined as

[TABLE]

with

[TABLE]

where $\Gamma_{1r}=\Gamma_{1r}(\boldsymbol{X},T)\in\mathcal{S}(\mathbb{R}^{2})~{}(\forall~{}T\geq 0)$, and $\Psi_{1r}=\Psi_{1r}(\boldsymbol{x})\in C^{\infty}(\mathbb{R}^{2})\cap L^{2}_{\mathbf{K}}(\mathbb{R}^{2}/\Lambda),~{}r\in\{1,\cdots,6\}$ have been defined in (5.7). By Poisson-Summation stated in (5.1), we can directly have

[TABLE]

Let us start to deal with ${G_{1}}_{\mathrm{D,I}}(\boldsymbol{x},t)$ at first. For any $\frac{1}{2}<\nu_{1}<1$,

[TABLE]

Then, by the formula (5.3), we denote that

[TABLE]

According to the statement (5.4), when $\mathbf{m}\neq(0,0)$, $|\mathbf{k}-\mathbf{K}|<\varepsilon^{\nu_{1}}$ as $\varepsilon>0$ sufficiently small, there exists a positive number $C>0$ such that

[TABLE]

Due to the conclusion (5.2), it follows that for any $0\leq\tau\leq\rho\varepsilon^{{}^{-1}}$,

[TABLE]

Since $\Psi_{1r}(\boldsymbol{x}),~{}r=1,~{}\cdots,~{}6$ are smooth and $\{\Phi_{b}(\boldsymbol{x};\mathbf{k})\}_{b\geq 1}$ is also complete orthogonal in $L^{2}(\Omega)$ for any fixed $\mathbf{k}\in\Omega^{*}$, we immediately obtain

[TABLE]

and choose $N>2$ to guarantee the above convergence for any $\tau\in[0,~{}\rho\varepsilon^{{}^{-1}}]$.

To estimate $\text{I}_{1}(\mathbf{k},\tau)$, we will utilize the expansion of $\Phi_{\pm}(\boldsymbol{x};\mathbf{k})$ near $\mathbf{k}=\mathbf{K}$ which is given in Proposition 2.3, i.e., for $|\mathbf{k}-\mathbf{K}|<\varepsilon^{\nu_{1}}\leq q_{0}$,

[TABLE]

Substituting (5.25) into $\text{I}_{1}(\mathbf{k},\tau)$ and by the fact all $\Gamma_{1r}(\cdot,\varepsilon\tau)\in\mathcal{S}(\mathbb{R}^{2})$ uniformly for $0\leq\tau\leq\rho\varepsilon^{{}^{-1}}$, we directly obtain

[TABLE]

In the above calculation, the key is to derive \big{\langle}\Phi_{j}(\boldsymbol{x}),~{}\Psi_{1r}(\boldsymbol{x})\big{\rangle}_{L^{2}(\Omega)},~{}j=1,~{}2,~{}r=1,\cdots,6. By the conclusions stated in Proposition 2.4, all non-vanishing leading order terms in (5.3) are listed

[TABLE]

Since $\boldsymbol{\alpha}(\varepsilon\boldsymbol{x},\varepsilon t)$ is the solution of Dirac equation (1.7), we have the following identity by Fourier transform,

[TABLE]

Then it gives that for any $\tau\in[0,~{}\rho\varepsilon^{{}^{-1}}]$ and $|\mathbf{k}-\mathbf{K}|<\varepsilon^{\nu_{1}}$,

[TABLE]

Thanks to (5.24) and (5.28), for all $t\in[0,~{}\rho\varepsilon^{{}^{-1}}]$, we have

[TABLE]

In the next, we show the estimate of $\|{G_{1}}_{\mathrm{D,II}}(\boldsymbol{x},t)\|_{L^{2}(\mathbb{R}^{2})}^{2}$. Recall that definition

[TABLE]

Due to the result stated in (5.5), if $|\mathbf{k}-\mathbf{K}|\geq\varepsilon^{\nu_{1}}$, there exists a constant $C$ independent of $\varepsilon$ such that for all $\mathbf{m}\in\mathbb{Z}^{2}$,

[TABLE]

Then by invoking (5.2) and Poisson-Summation formula (5.3), we can obtain

[TABLE]

Thus, to ensure the convergence of (5.3) and (5.31), we need to choose $\frac{1}{2}<\nu_{1}<1$ i.e., $\nu_{1}=1-\frac{\nu}{2}$ which indicates $0<\nu<1$ should not be zero, and let $N$ be large enough such that for any $t\in[0,~{}\rho\varepsilon^{{}^{-1}}]$,

[TABLE]

The above (5.32) and (5.33) imply

[TABLE]

In order to estimate ${G_{1}}_{\mathrm{D^{C},I}}(\boldsymbol{x},t)$, let us fist introduce that

[TABLE]

where $q_{1}$ is chosen in Proposition 2.3 such that if $|\mathbf{k}-\mathbf{K}|<q_{1}$ and $b\notin\{+,-\}$, the following uniform bound holds

[TABLE]

and in addition, $\{E_{b}(\mathbf{k})\}_{b\geq 1,~{}|\mathbf{k}-\mathbf{K}|<q_{1}}$ have a positive lower bound.

By carrying out the integration by part to $\tilde{G_{1}}_{b}(\mathbf{k},t)$ in $t$, it yields that

[TABLE]

and similarly,

[TABLE]

Utilizing (5.36) and the unform lower bound of $E_{b}(\mathbf{k})$ when $|\mathbf{k}-\mathbf{K}|<q_{1}$, we immediately conclude the following estimate holds

[TABLE]

Then by Minkowski’s integral inequality, for any $t\in[0,~{}\rho\varepsilon^{{}^{-1}}]$ we have

[TABLE]

We need a new approach to estimate ${G_{1}}_{\mathrm{D^{C},II}}(\boldsymbol{x},t)$. Indeed, the technique of integration by part in the previous estimate will be invalid, since the spectral band $E_{b}(\mathbf{k})$ varies unclearly when $|\mathbf{k}-\mathbf{K}|\geq q_{1}$ and $|\mathbf{k}|\geq\delta_{b,1}\frac{|\mathbf{K}|}{2}$ in $\Omega^{*}$. Here the strategy is similar to that of verifying ${G_{1}}_{\mathrm{D,II}}$. Main difference is that we need to control under the $H^{s}$ norm, thus a sharper estimate is carried out by applying $(I+\mathcal{L}^{A})$ on $\Phi_{b}(\boldsymbol{x};\mathbf{k})$ as many times as required. Namely, for any integer $M\geq 0$,

[TABLE]

Note that (I+\mathcal{L}^{A})^{M}\big{(}\Gamma_{1r}(\varepsilon\boldsymbol{x},\varepsilon\tau)\cdot\Psi_{1r}(\boldsymbol{x})\big{)} is a summation of terms in the form $\Upsilon(\varepsilon\boldsymbol{x},\varepsilon\tau)\Theta(\boldsymbol{x})$ where $\Upsilon(\cdot,\varepsilon\tau)\in\mathcal{S}(\mathbb{R}^{2})$ includes derivatives of $\Gamma_{1r}(\varepsilon\boldsymbol{x},\varepsilon\tau)$, and $\Theta(\boldsymbol{x})\in C^{\infty}(\mathbb{R}^{2})\cap L^{2}_{\mathbf{K}}(\mathbb{R}^{2}/\Lambda)$ contains derivatives of $A(\boldsymbol{x})$ and $\Psi_{1r}(\boldsymbol{x})$. According to (5.3) in Proposition 5.1, if $\mathbf{k}\in\Omega^{*}$ and $|\mathbf{k}-\mathbf{K}|\geq q_{1}$, there exists a constant $C>0$ such that

[TABLE]

Therefore, by applying (5.1) and (5.2), we can conclude that for any $\mathbf{k}\in\Omega^{*},~{}|\mathbf{k}-\mathbf{K}|\geq q_{1}$ and $\tau\in[0,~{}\rho\varepsilon^{{}^{-1}}]$,

[TABLE]

Thus, one can figure out for any $t\in[0,~{}\rho\varepsilon^{{}^{-1}}]$

[TABLE]

Here we require $N>2$ to ensure the convergence of the double summation. This leads to the last estimation as follows

[TABLE]

where we choose $M>(s+1)/2$.

Finally, (5.34), (5.3), (5.3) and Lemma 5.2 imply that for any $s\geq 0$, $\rho>0$, $0<\nu<1$ and $\varepsilon>0$ sufficiently small

[TABLE]

This completes the proof of Proposition 4.3 and therewith Theorem 4.1.

$\square$

Acknowledgments: The authors would thank Prof. Michael I. Weinstein for proposing this interesting problem and for his useful suggestions. This work was supported by National Natural Science Foundation of China NSFC grants $\#11871299$.

A Appendix: Wave packets with data spectrally localized near Dirac points

In this appendix, we consider the case of $\kappa(\varepsilon\boldsymbol{x})\equiv 0$ and aim to obtain an asymptotic solution that would be valid over the time scale up to the order $\mathcal{O}(\frac{1}{\varepsilon^{{}^{2^{-}}}})$, which is parallel to the main result in [20]. To put the analysis in the context, consider the following wave packet equations with data spectrally localized near Dirac points

[TABLE]

with the same initial conditions (1.5).

We obtain a similar asymptotic solution in the form of

[TABLE]

Here the envelopes $\alpha_{j}(\varepsilon\boldsymbol{x},\varepsilon t),~{}j=1,~{}2$ satisfy the massless Dirac equation

[TABLE]

with initial condition $\alpha_{1}(\boldsymbol{X},0)=\alpha_{10}(\boldsymbol{X}),~{}\alpha_{2}(\boldsymbol{X},0)=\alpha_{20}(\boldsymbol{X})$.

To conclude, we can obtain the following result,

Theorem A.1**.**

Supposed that $A(\boldsymbol{x})$ is a honeycomb structured material weight defined in Definition 2.1, $\Phi_{1}(\boldsymbol{x})$, $\Phi_{2}(\boldsymbol{x})$ are the eigenfunctions associated with the Dirac point $(\mathbf{K},E_{D})$ given in Definition 2.2, $\alpha_{10}(\boldsymbol{X})$, $\alpha_{20}(\boldsymbol{X})$ are Schwartz functions of $\boldsymbol{X}\in\mathbb{R}^{2}$. Then problem (A.1) (1.5) has a unique solution of the form (A.2), where $\alpha_{j}(\boldsymbol{X},T),~{}j=1,~{}2$ are the solution to the massless system (A.3), and for any $s\geq 0$, $\rho>0$, $0<\nu<1$, a sufficiently small positive parameter $\varepsilon$,

[TABLE]

Here $C$ does not depend on $\varepsilon$.

Firstly, the well-posedness of the massless Dirac equation (A.3) for $T\in[0,\infty)$ can be obtained by the Fourier transform, see [20]. Since the error estimate is quite similar to that in the proof of Theorem 4.1, we omit the detailed repeated calculations by just sketching out the main idea. After substituting (A.2) into equation (A.1) and decomposing $\eta(\boldsymbol{x},t)$ into its Floquet-Bloch components, we obtain a similar equation (4) (4) without the implicit term $\tilde{\mathcal{Q}}_{b}[\eta]$. By carefully dealing with (4.16)-(4.18), we can improve the estimates as follows

[TABLE]

Thus Theorem A.1 is asserted.

B Appendix: Proof of Proposition 3.1

The proof of Proposition 3.1 is given in this appendix. We first list the standard result on the well-posedness of linear symmetric hyperbolic system, see for instance [26, 35].

Let $U=U(\boldsymbol{x},t)$ be the unknown $n$-vector-valued function from $(\boldsymbol{x},t)\in\mathbb{R}^{d}\times[0,+\infty)$ to $\mathbb{R}^{n}$, and $A_{1}(\boldsymbol{x},t),~{}\cdots,~{}A_{d}(\boldsymbol{x},t),~{}B(\boldsymbol{x},t)$ are $n\times n$-matrix-valued functions. Consider the first-order partial differential equations of the form

[TABLE]

with initial condition $U(\boldsymbol{x},0)=U_{0}(\boldsymbol{x})$. The global existence of the solution is shown in the following Proposition.

Proposition B.1**.**

Let $s\in\mathbb{N}$ and $s>\frac{d}{2}+1$. Assume that $A_{1},\cdots,A_{d}\in C_{b}^{s+1}[~{}or~{}C_{b}^{\infty}]$ are Hermitian matrices, $B\in C_{b}^{s+1}[~{}or~{}C_{b}^{\infty}]$ and $U_{0}\in H^{s}[~{}or~{}\cap C^{\infty}]$. Then there exists a unique solution U\in C^{1}\big{(}\mathbb{R}^{d}\times[0,+\infty)\big{)}[~{}or~{}C^{\infty}\big{(}\mathbb{R}^{d}\times[0,+\infty)\big{)}], and for any $0<\rho<\infty$,

[TABLE]

Moreover, there exist $C_{1},~{}C_{2}>0$ such that

[TABLE]

To prove Proposition 3.1, we rewrite the Dirac equation (1.7) in the compact form

[TABLE]

where the coefficients $A_{1}=-\frac{v_{{}_{F}}}{2\sqrt{E_{D}}}\sigma_{1},~{}A_{2}=\frac{v_{{}_{F}}}{2\sqrt{E_{D}}}\sigma_{2}$ are Hermitian matrices and $B(\boldsymbol{X})=-i\frac{\vartheta\kappa(\boldsymbol{X})}{2\sqrt{E_{D}}}\sigma_{3}$ is a smooth bounded matrix-valued function.

Note that the initial value \boldsymbol{\alpha}_{0}(\boldsymbol{X})=\big{(}\alpha_{10}(\boldsymbol{X}),\alpha_{20}(\boldsymbol{X})\big{)}^{T} is in Schwartz space. The first conclusion (3.1) of Proposition 3.1 is just a direct consequence of Proposition B.1. However, it requires a more delicate estimate to derive the second result (3.2) in Proposition 3.1.

For any integer $N\geq 0$ and the multi-indices $\mathbf{n}\in\mathbb{N}^{2}$ satisfying $|\mathbf{n}|\leq N$, define

[TABLE]

According to the Dirac equation (1.7), one can deduce $V_{N}$ satisfies

[TABLE]

with the initial value

[TABLE]

where

[TABLE]

are both block diagonal matrices, and $\mathscr{B}(\boldsymbol{X})$ is a lower triangular matrix with each block entry being a linear combination of $\partial^{\mathbf{m}}_{\boldsymbol{X}}B(\boldsymbol{X}),~{}\mathbf{m}\in\mathbb{N}^{2},~{}|\mathbf{m}|\leq N$, and thus smooth and bounded.

Therefore, for any $M\in\mathbb{N}$, $(1+|\boldsymbol{X}|^{2})^{\frac{M}{2}}V_{N}$ satisfies the following system

[TABLE]

with initial values

[TABLE]

where $\widetilde{\mathscr{B}}(\boldsymbol{X})=\mathscr{B}(\boldsymbol{X})-\mathscr{A}_{1}\frac{MX_{1}}{1+|\boldsymbol{X}|^{2}}-\mathscr{A}_{2}\frac{MX_{2}}{1+|\boldsymbol{X}|^{2}}$ is also smooth and bounded.

By Proposition B.1, we can conclude for any $s>2,~{}0<\rho<\infty$,

[TABLE]

and further for any $T\geq 0$,

[TABLE]

Then (B.11) implies $(1+|\boldsymbol{X}|^{2})^{\frac{M}{2}}V_{N}$ is in H$\ddot{o}$lder space by Sobolev embedding theorem, and for any $\rho>0$,

[TABLE]

If $l\in\mathbb{N}$, one can observe by induction, $\partial_{T}^{l}\boldsymbol{\alpha}$ is a linear combination of $\partial_{\boldsymbol{X}}^{\mathbf{m}}\boldsymbol{\alpha},~{}0\leq|\mathbf{m}|\leq l$, with the smooth bounded coefficients. Specifically, $\partial_{T}^{2}\alpha_{1},~{}\partial_{T}^{2}\alpha_{2}$ satisfy an uncoupled form:

[TABLE]

In conclusion, one can get for any $l,~{}M\in\mathbb{N},~{}\mathbf{n}\in\mathbb{N}^{2}$, there exists a positive constant $C$ such that

[TABLE]

This completes the proof.

$\square$

References