Wave-packet dynamics in multilayer phosphorene

We investigate the dynamics of Gaussian wave packets in multilayer black phosphorus (BP). Time-dependent average position and velocity are calculated analytically and numerically by using a continuum model and a method based on the split-operator technique, respectively. By analyzing the wave-packet trajectories with nonvanishing initial momentum along $y$ direction, we observed transient spatial oscillations due to the effect known as zitterbewegung (ZBW). We also demonstrated, based on the Heisenberg picture and by the calculation of the velocity operators, that the trembling motion along the $y$ direction at small times is unavoidable even for null initial momentum. We verified that the ZBW is directly related to the splitting of the wave packet into two parts moving with opposite velocities, as similar to graphene, and the linear dependence on $k_y$ in out-of-diagonal terms in the Hamiltonian. In addition for phosphorene systems, the two portions of the propagated wave packets have an asymmetric shape for unbalanced $\left({[1,0]}^T\right)$ and phased different $\left({[1,i]}^T\right)$ initial pseudospin components, playing an important role in the amplitude, frequency and duration time of the transient oscillations. Electrons in the vicinity of the Fermi energy traveling in $N$-layer phosphorene also exhibit qualitatively similar trembling motion and transient character as found in the monolayer case, except by the oscillation phase difference and final group velocity achieved after the transient behavior. As a consequence of the anisotropy on the $N$-layer BP energy bands, effective masses and group velocities along the $x$ and $y$ directions, the wave packet propagates nonuniformly along the different directions and deforms into an elliptical shape. By comparing our analytical results with those ones obtained by the split-operator technique, we verified a good qualitative and quantitative agreement between them, except for very larger values of wave vector and after long time steps.

Figure 1

(a) Top view of lattice structure of AB-stacked $N$-layer BP system, emphasizing the orientations of the lattice adopted in this work and the four sublattices: $A$ and $B$ at bottom sublayer (purple symbols), and $C$ and $D$ at top sublayer (brown symbols). The $x$ and $y$ coordinates correspond to the zigzag and armchair directions, respectively, and $z$ direction is the out-of-plane direction. (b) (Left) Lowest electronic energy band obtained by diagonalizing the Hamiltonian (1) with $n=N$ for monolayer (black solid curve), bilayer (red dashed curve), trilayer (blue dot curve), and tetralayer (green short-dash curve) phosphorene. (Right) Initial wave-packet energy for the corresponding initial wave vector assumed here in the wave-packet simulation. ${\theta }_k$ as a function of (c) the polar angle $\beta$ for fixed momentum vector and (d) the momentum vector for fixed polar angle $\beta$. (e) Momentum value for ${\theta }_{\max }$ as a function of the polar angle $\beta$, i.e., $\left|k\right|$ in which ${\theta }_k$ has a maximum value as emphasized in the orange dashed line in panel (d).

Figure 2

[(a) and (b)] Average position and [(c) and (d)] expectation value of the velocity for trajectories of (a) and (b), respectively, of a Gaussian wave packet of width $d=100$ Å as a function of time for the case $c_1=1$ and $c_2=0$. [(a) and (c)] Wave packet propagating in a monolayer phosphorene sheet $\left(N=1\right)$ with different initial central wave vectors: $k_0=0.01$ (black), 0.05 (red) and $0.1{\AA }^{-1}$ (blue). [(b) and (d)] Wave packet propagating in multilayer phosphorene $(N=1$, 2, 3, 4) for fixed wave vector $k_0=0.1{\AA }^{-1}$. The wave packet starts at $(x_0,y_0)=(0,0)$ Å. The inset shows an enlargement of the physical averages [(a) and (b)] for first time steps in order to emphasize the different oscillation amplitudes and oscillation frequencies, and (d) nearby the time steps in which the velocities achieve constant values.

Figure 3

Time evolution of the electronic wave packet for the case ${[c_1,c_2]}^T={[1,0]}^T$ with $\left|k\right|=0.05{\AA }^{-1}$, corresponding the average position shown by the red curve in Fig. 2. Snapshots at (a) $t=20$, (b) 30, (c) 40 and (d) 50 fs.

Figure 4

The same as in Fig. 2, but now for the case $c_1=1$ and $c_2=1$.

Figure 5

The same as in Fig. 3, but now for the case $c_1=1$ and $c_2=1$.

Figure 6

The same as in Fig. 2, but now for the case $c_1=1$ and $c_2=i$.

Figure 7

The same as in Fig. 3, but now for the case $c_1=1$ and $c_2=i$.

Figure 8

Comparison between the results of the average position of a Gaussian wave packet of width $d=100$ Å as a function of time obtained by (solid curves) the split-operator technique derived in Sec. 2c and (dashed curves) analytical calculations derived in Sec. 2b: [(a) and (b)] the case $c_1=1$ and $c_2=1$ and [(c) and (d)] the case $c_1=1$ and $c_2=i$. [(a) and (c)] Wave packet propagating in a monolayer phosphorene sheet $\left(N=1\right)$ with different wave vectors: $k_0=0.01$ (black), 0.05 (red), and $0.1{\AA }^{-1}$ (blue). [(b) and (d)] Wave packet propagating in multilayer phosphorene $(N=1$, 2, 3, 4) for fixed wave vector $k_0=0.1{\AA }^{-1}$.

Figure 9

(a) Average position and [(b) and (c)] expectation value of the velocity along $y$ direction as a function of time for the case ${[c_1,c_2]}^T={[1,0]}^T$ (black), ${[c_1,c_2]}^T={[1,1]}^T$ (red), and ${[c_1,c_2]}^T={[1,i]}^T$ (blue). (c) A magnification of the results in (b) for large time steps showing the oscillatory behavior of $v_y$. (Right) Contour plots of the squared modulus of the wave function at (I) $t=390$ and (II) 391 fs, and a zoom emphasized by the dashed curves showing the isosurfaces at the two time steps.