# Bursty Star Formation Naturally Explains the Abundance of Bright Galaxies at
Cosmic Dawn

Guochao Sun  CIERA and Department of Physics and Astronomy, Northwestern
University, 1800 Sherman Ave, Evanston, IL 60201, USA  Claude-André Faucher-
Giguère  CIERA and Department of Physics and Astronomy, Northwestern
University, 1800 Sherman Ave, Evanston, IL 60201, USA  Christopher C. Hayward
Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue,
New York, NY 10010, USA  Xuejian Shen  TAPIR, California Institute of
Technology, Pasadena, CA 91125, USA  Department of Physics & Kavli Institute
for Astrophysics and Space Research, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Andrew Wetzel  Department of Physics & Astronomy,
University of California, Davis, CA 95616, USA  Rachel K. Cochrane  Center for
Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY
10010, USA

###### Abstract

Recent discoveries of a significant population of bright galaxies at cosmic
dawn $\left(z\gtrsim 10\right)$ have enabled critical tests of cosmological
galaxy formation models. In particular, the bright end of the galaxy UV
luminosity function (UVLF) appears higher than predicted by many models. Using
approximately 25,000 galaxy snapshots at $8\leq z\leq 12$ in a suite of FIRE-2
cosmological “zoom-in” simulations from the Feedback in Realistic Environments
(FIRE) project, we show that the observed abundance of UV-bright galaxies at
cosmic dawn is reproduced in these simulations with a multi-channel
implementation of standard stellar feedback processes, without any fine-
tuning. Notably, we find no need to invoke previously suggested modifications
such as a non-standard cosmology, a top-heavy stellar initial mass function,
or a strongly enhanced star formation efficiency. We contrast the UVLFs
predicted by bursty star formation in these original simulations to those
derived from star formation histories (SFHs) smoothed over prescribed
timescales (e.g., 100 Myr). The comparison demonstrates that the strongly
time-variable SFHs predicted by the FIRE simulations play a key role in
correctly reproducing the observed, bright-end UVLFs at cosmic dawn: the
bursty SFHs induce order-or-magnitude changes in the abundance of UV-bright
($M_{\mathrm{UV}}\lesssim-20$) galaxies at $z\gtrsim 10$. The predicted
bright-end UVLFs are consistent with both the spectroscopically confirmed
population and the photometrically selected candidates. We also find good
agreement between the predicted and observationally inferred integrated UV
luminosity densities, which evolve more weakly with redshift in FIRE than
suggested by some other models.

galaxies: formation – galaxies: evolution – galaxies: star formation –
galaxies: high-redshift

††journal: ApJL††software: BPASS (Eldridge et al., 2017), GizmoAnalysis
(Wetzel et al., 2016; Wetzel & Garrison-Kimmel, 2020), hmf (Murray et al.,
2013)

##  1 Introduction

For the first time, the James Webb Space Telescope (JWST) has unlocked the
door to a population-level analysis of galaxies well into the era of cosmic
dawn (for a review of key high-redshift science themes of JWST, see Robertson,
2022). Following its discovery of an unexpectedly high abundance of UV-bright,
massive galaxy candidates at redshift $z\gtrsim 10$ (e.g., Finkelstein et al.,
2022; Naidu et al., 2022a; Donnan et al., 2023; Harikane et al., 2023b; Yan et
al., 2023), there is a long list of intriguing questions to be answered about
how to interpret these observations. What is the true nature (redshift, mass,
metallicity, age, etc.) of these bright galaxies? If they are truly massive
galaxies at cosmic dawn, what makes it possible for them to have formed so
early? Are these observations in significant tension with the standard
$\Lambda$CDM cosmological model? Observational and theoretical investigations
into these questions are being actively pursued in a large body of recent
literature from different perspectives, including the purity of high-$z$
candidates (Naidu et al., 2022b; Arrabal Haro et al., 2023; Curtis-Lake et
al., 2023; Furlanetto & Mirocha, 2023; Zavala et al., 2023), the physics of
star formation in high-$z$ galaxies (Dekel et al., 2023; Mirocha & Furlanetto,
2023; Robertson et al., 2023; Qin et al., 2023; Sipple & Lidz, 2023; Trinca et
al., 2023), the implications of high-$z$ observations for the cosmological
model (Boylan-Kolchin, 2023; Hassan et al., 2023; Melia, 2023), and so forth.

While spectroscopic follow-up studies for many of the galaxy candidates are
still ongoing, conservative lower limits on the bright end of UV luminosity
function (UVLF) and the integrated UV luminosity density at $z\gtrsim 10$
derived from the existing, spectroscopically confirmed samples have already
suggested milder redshift evolution towards $z>10$ than expected by many
theoretical models (e.g., Harikane et al., 2023a). Such a higher-than-expected
abundance of bright galaxies based on secure redshifts is consistent with
earlier studies based on photometrically selected samples, thus calling for a
re-examination of the theoretical landscape of galaxy formation at cosmic
dawn111Some recent studies found that galaxies with properties similar to
observed ones could be reproduced in simulations (e.g., Keller et al., 2023;
McCaffrey et al., 2023). However, these studies did not directly model the
UVLF and compare it with available JWST measurements.. Several physical
mechanisms have been considered to explain a high abundance of bright galaxies
at high redshifts. For example, a higher star formation efficiency (SFE)
resulting from less efficient feedback regulation could boost the UV-bright
galaxy abundance by forming more stars per unit baryon (Dekel et al., 2023;
Harikane et al., 2023b), whereas a more top-heavy initial mass function (IMF)
of the stellar population could similarly lead to more bright galaxies by
creating more UV photons per unit stellar mass formed (Inayoshi et al., 2022;
Yung et al., 2023). A conspiracy between the redshift evolution of dust
attenuation and the abundance of massive halos at high $z$ could also
potentially allow the bright-end UVLF and UV luminosity density to evolve
relatively mildly (Ferrara et al., 2023; Mirocha & Furlanetto, 2023), although
such a coincidence would not by itself explain the correct absolute abundance
of bright galaxies. A number of studies have also examined the possibility
that the high abundance of early massive galaxies implies physics beyond the
standard $\Lambda$CDM cosmology, such as a modified primordial power spectrum
(Hirano & Yoshida 2023; Padmanabhan & Loeb 2023; Parashari & Laha 2023; though
see Sabti et al. 2023), primordial non-Gaussianity (Biagetti et al., 2023), or
alternative dark matter models (Bird et al., 2023; Dayal & Giri, 2023; Gong et
al., 2023).

Another promising avenue to elevate the abundance of bright galaxies is the
strong time variability (“burstiness”) of star formation. In recent years,
several different galaxy formation simulations have predicted that the star
formation rate (SFR) is highly time-variable in low-mass galaxies (e.g.,
Hopkins et al., 2014; Domínguez et al., 2015; Muratov et al., 2015; Sparre et
al., 2017; Pallottini & Ferrara, 2023). The prediction of bursty star
formation appears generic to codes that resolve the clustering of supernovae
in the interstellar medium (Hu et al., 2023). The simulations predict that
bursty star formation is especially common in low-mass galaxies, likely due to
the shallow potential wells which allow clumpy, cold inflows and outflows to
drive repeated inflow-star formation-outflow cycles (Stern et al., 2021;
Gurvich et al., 2023; Byrne et al., 2023; Hopkins et al., 2023). Since low-
mass galaxies dominate at high redshift, we expect the implications of bursty
star formation on the UVLF to be particularly important in this regime (e.g.,
Furlanetto & Mirocha, 2022). Indeed, evidence for an increased level of bursty
star formation has emerged from recent JWST observations of cosmic dawn
galaxies (e.g., Dressler et al., 2023; Endsley et al., 2023; Looser et al.,
2023a, b). As pointed out in recent theoretical studies (Mason et al., 2023;
Mirocha & Furlanetto, 2023; Shen et al., 2023; Muñoz et al., 2023), an
increased level of UV variability sourced by bursty star formation can give
rise to more UV-bright galaxies due to the Eddington bias, which flattens the
bright end of the UVLF. In this case, the observed UVLFs at $z\gtrsim 10$
could potentially be explained by bursty star formation combined with “normal”
SFE and production efficiency of UV photons. While bursty star formation can
in principle enhance the abundance of bright galaxies, it remains to be shown
whether the enhancement is sufficient to reproduce the observed bright-end of
the UVLF in a self-consistent galaxy formation model, such as those provided
by hydrodynamic simulations.

In this Letter, we use a suite of cosmological “zoom-in” simulations from the
Feedback in Realistic Environments (FIRE) project222See the FIRE project
website: http://fire.northwestern.edu. to investigate the effects of bursty
star formation on the UVLF at $8\leq z\leq 12$. In these simulations, the SFR
variability arises self-consistently from the modeling of standard stellar
feedback processes. It is noteworthy that these simulations — generated before
the launch of JWST — were in particular not in any way tuned to match recent
observations. Moreover, the simulations use exactly the same FIRE-2 code
(Hopkins et al., 2018) that has been used to evolve large sets of simulated
galaxies all the way to $z=0$ and demonstrated to produce broadly realistic
galaxy properties down to the present time (e.g. Wetzel et al., 2023, and
references therein). This is in contrast with many other simulations of cosmic
dawn galaxies, in which the simulations are stopped at high redshift and for
which we therefore do not know how the feedback model performs at lower
redshifts. We show that the FIRE-2 simulations produce an excellent match to
the UVLF recently measured by JWST during cosmic dawn, and that the time
variability of star formation plays an important role in explaining the
observations at the bright end. These results constitute an important test of
the feedback model and highlight the importance of considering the variability
of star formation when modeling high-$z$ observations.

Throughout the Letter, we adopt a flat $\Lambda$CDM cosmology consistent with
Planck Collaboration et al. (2020), and all magnitudes are quoted in the AB
system (Oke & Gunn, 1983).

##  2 Simulations and Analysis Methods

###  2.1 The Simulations

In this Letter, we analyze the same set of simulations as recently studied by
Sun et al. (2023a), which is a subset of the High-Redshift suite (Ma et al.,
2018a, b, 2019) of the FIRE-2 cosmological zoom-in simulations (Hopkins et
al., 2018). The FIRE-2 simulations use the GIZMO code with its meshless-finite
mass (MFM) hydro solver (Hopkins, 2015), and include multiple channels of
stellar feedback to regulate star formation. Star formation occurs in dense
molecular gas ($n_{\mathrm{H}}>1000\,\mathrm{cm^{3}}$) that is self-
gravitating and self-shielding. The stellar feedback mechanisms implemented
include: (1) energy, momentum, mass, and metal injection from core collapse
and Type Ia supernovae and winds from OB and AGB stars, (2) photoionization
and photoelectric heating, and (3) radiation pressure. A redshift-dependent
but homogeneous ionizing background is also included following Faucher-Giguère
et al. (2009).333The version of the ionizing background used in these
simulations reionizes the universe at $z_{\rm reion}\approx 10$, which is
earlier than the mid-point of reionization of $z_{\rm reion}\approx 8$ favored
by more recent observational constraints (e.g., Planck Collaboration et al.,
2020; Faucher-Giguère, 2020). However, our main results focus on the bright
end of the UVLF, which arises from relatively massive halos, whereas the
suppression of galaxy formation due to heating by the ionizing background
primarily affects low-mass halos ($M_{\rm h}\lesssim 10^{9}$ $M_{\odot}$; e.g.
Gnedin, 2000; Noh & McQuinn, 2014). Moreover, an earlier reionization redshift
implies that in the present simulations, galaxy formation is suppressed
starting earlier in the small halos, so adopting a more up-to-date
reionization model would (if anything) enhance the predicted UV luminosity
density. Similar arguments apply to other IGM heating processes. The baryonic
(dark matter) mass resolution of the set of simulations considered in this
work is $m_{\mathrm{b}}=7\times 10^{3}\,M_{\odot}$ ($m_{\rm DM}=4\times
10^{4}\,M_{\odot}$), except for the simulations z5m11a and z5m11b, which have
$m_{\mathrm{b}}\approx 1\times 10^{3}\,M_{\odot}$ ($m_{\rm DM}=5\times
10^{3}\,M_{\odot}$). The gravitational softenings are fixed in physical units
to $\epsilon_{\mathrm{DM}}=42$ pc for the dark matter and
$\epsilon_{\mathrm{star}}=2.1$ pc for stars. The gravitational softenings are
adaptive for gas, with a minimum of $\epsilon_{\mathrm{b}}=0.42$ pc. This is,
again, with the exception of z5m11a and z5m11b (see Figure 4 for a list of
simulation IDs considered in this work), which have
$\epsilon_{\mathrm{DM}}=21$ pc, $\epsilon_{\mathrm{star}}=1.4$ pc, and
$\epsilon_{\mathrm{b}}=0.28$ pc.

Part of the High-Redshift suite of simulations was presented and analyzed in
detail by Ma et al. (2018a, b) for the predicted properties of the simulated
galaxy population at $5\leq z\leq 12$, including sizes, morphologies, scaling
relations, and number statistics measured by the stellar mass and luminosity
functions. In this follow-up analysis of Ma et al. (2018b) motivated by recent
JWST observations of the abundance of galaxies at $z\gtrsim 10$, we follow
closely the methodology adopted in Ma et al. (2018b) for fair comparisons, but
the sample size of high-$z$, massive galaxies has been substantially increased
to better determine the bright-end behavior of galaxy UVLFs at cosmic dawn.
Below, we will only briefly summarize the key information about the sample of
simulated galaxies pertinent to the analysis presented here. We refer
interested readers to the aforementioned papers for further details about the
FIRE-2 simulations and the High-Redshift suite.

For a robust analysis of UVLFs at their bright end, we build a maximum
possible sample size of massive galaxies by making use of all the zoom-in
simulations available at each redshift above the ending redshift
$z_{\mathrm{end}}$. In each zoom-in region, we consider all the well-resolved
halos444Halos containing at least $10^{4}$ particles in total and
uncontaminated by low-resolution particles are considered “well-resolved”.
that host a central galaxy, rather than the one hosting just the most massive,
primary galaxy (typically near the center of the zoom-in region). Following Ma
et al. (2018b), we define galaxies based on catalogs of halos identified with
the Amiga Halo Finder (AHF; Knollmann & Knebe 2009). The radius
$R_{\mathrm{max}}$ at which the halo rotation curve reaches maximum is used to
define a galaxy by incorporating star particles within $R_{\mathrm{max}}/3$
and excluding the contamination from subhalos outside $R_{\mathrm{max}}/5$. We
restrict the scope of our UVLF analysis to halos with mass
$M_{\mathrm{h}}>10^{7.5}\,M_{\odot}$ in snapshots at $8\leq z\leq 12$ because
most of the recent UVLF measurements at $z<8$ with JWST can be well explained
by previous theoretical predictions and a sufficiently constraining sample of
spectroscopically-confirmed galaxies is not available at $z>12$ (Harikane et
al., 2023a). In Appendix A, we illustrate how the halo/galaxy sample is
constructed with (snapshots of) the 26 individual zoom-in simulations, which
build up a total sample of $\approx 25,000$ galaxy snapshots over $8\leq z\leq
12$. For all simulations, snapshots are saved at a cadence of every 10–20 Myr.

###  2.2 Processing of the Simulations

We process the simulated galaxy sample in order to arrive at their 1600 Å UV
magnitudes $M_{\mathrm{UV}}$, following Sun et al. (2023a). Templates of
binary, single-stellar-population (SSP) spectra from BPASS v2.1 (Eldridge et
al., 2017) are interpolated and applied to star particles according to their
stellar age and metallicity, assuming a Kroupa IMF (Kroupa, 2001). Including
nebular (continuum) emission can in principle augment both the UV emissivity
and variability (Byler et al., 2017), although we opt to ignore it here as
nebular emission is not expected to strongly affect the measurement of
$M_{\mathrm{UV}}$, especially when compared with effects of SFR variations.
Two notable differences from Sun et al. (2023a) exist, though, for the
treatment of (1) the connection between $M_{\mathrm{UV}}$ and the SFH and (2)
dust attenuation, on which we elaborate below.

####  2.2.1 Bursty vs Smoothed Star Formation Histories

At cosmic dawn, an increased SFR variability can strongly modulate the
observed number statistics of galaxies. To assess the impact of bursty star
formation on the $M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation and thus the
UVLF, we consider two contrasting scenarios to model $M_{\mathrm{UV}}$.

The baseline scenario, which we refer to as “bursty”, assumes that the SFH of
each galaxy in our sample is exactly as predicted by the simulations and thus
$M_{\mathrm{UV}}$ can be derived by summing up the spectral emissivities of
all star particles of the galaxy at a given snapshot according to their age
and metallicity, as in Sun et al. (2023a). This is the approach most faithful
to the SFHs predicted by the simulations. In this approach, $M_{\mathrm{UV}}$
naturally inherits the burstiness predicted by the simulations — as the SFR
varies, the UV 1600 Å luminosity of the galaxy also fluctuates accordingly
because most of FUV continuum emission is sourced by the massive, short-lived
stars formed. As a result, a bursty SFH imprints significant stochasticity in
$M_{\mathrm{UV}}$ at a fixed stellar or halo mass.

In the contrasting scenario, which we refer to as “smoothed”, we artificially
reduce the impact of bursty SFH on the evaluation of $M_{\mathrm{UV}}$ by
redistributing the ages of star particles (while retaining their
metallicities). Specifically, we first define a smoothing kernel of duration
$\tau_{\mathrm{SF}}$ Myr and bin star particles using their star formation
times into time bins of width $\tau_{\mathrm{SF}}$. We then redistribute the
ages of the star particles in individual bins such that the stellar mass forms
at a nearly constant rate by enforcing evenly-distributed star formation times
within each bin. This redistribution of stellar ages effectively smooths the
SFH and reduces to the “bursty” case for a sufficiently small
$\tau_{\mathrm{SF}}$. Notably, unlike some previous work where effects of
varying the UV variability on UVLFs are studied assuming a fixed mean/median
$L_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation (e.g., Mirocha & Furlanetto, 2023;
Shen et al., 2023), our method by its nature conserves the total amount of
cosmic star formation such that the two scenarios differ only in terms of the
short-timescale SFR variability and its impact on the UV emissivity.

####  2.2.2 Dust Attenuation

Observations have shown compelling evidence of early chemical enrichment and
the production of non-negligible dust in galaxies at $z\gtrsim 7$ (Tamura et
al., 2019; Fudamoto et al., 2021; Witstok et al., 2023). A reasonable
treatment of dust attenuation is therefore needed for our predictions of the
UVLF at cosmic dawn, especially at the bright end because massive
(intrinsically UV-bright) galaxies generally contain more dust.

To estimate the effect of dust attenuation on $M_{\mathrm{UV}}$, we employ an
empirical model motivated by an up-to-date measurement of the
$\beta_{\mathrm{UV}}$–$M_{\mathrm{UV}}$ (color–magnitude) relation at $z>8$ by
Cullen et al. (2023) using a combination of JWST and ground-based
observations555See also Topping et al. (2023), who find slightly steeper
$\beta_{\mathrm{UV}}$ that steepens with increasing redshift from $z\sim
6$–12.. We combine the best-fit relation
$\beta_{\mathrm{UV}}=-0.17M_{\mathrm{UV}}+5.40$ with the attenuation–UV slope
relation, $A_{\mathrm{UV}}=0.48(\beta_{\mathrm{UV}}+2.62)$, determined from
$z\approx 5.5$ galaxies observed in the ALPINE survey (Fudamoto et al. 2020;
see also Reddy et al. 2018). While an extrapolation in redshift is involved,
this best-fit relation from ALMA observations represents a state-of-the-art
empirical baseline for estimating dust attenuation properties at cosmic dawn,
which should suffice for the purpose of this work. We neglect the scatter
around these mean relations given its small impact on $M_{\mathrm{UV}}$ and
caution that results with dust attenuation included that follow should be
taken as rough estimates only. The validity of these simplistic treatments can
be tested with simulations with detailed dust radiative transfer (Cochrane et
al., 2019, 2022; Ma et al., 2019; Vogelsberger et al., 2020; Shen et al.,
2022) and multi-wavelength observations (Akins et al., 2023; Bakx et al.,
2023), which are left for future work. We note, though, that at $z>10$ the
difference between UVLFs with and without dust attenuation is predicted to be
very small in our model (see Figure 2), such that uncertainties in the
treatment of dust should not affect our results significantly.

###  2.3 Estimating the UVLF from Zoom-in Simulations

Using UV magnitudes derived for the sample of simulated galaxies binned into
redshift bins of width $\Delta z=\pm 0.5$, we calculate the UVLF through a
convolution with the halo mass function (HMF) following the “HMF-weighting”
method introduced by Ma et al. (2018b). This method has been verified to
provide robust estimates of the UVLF from galaxy samples drawn from zoom-in
simulations, so we only summarize briefly here. First, in narrow halo mass and
redshift bins, we count the number of simulated halos $N_{\mathrm{S}}$ from
the sample and compute the expected number of halos $N_{\mathrm{E}}$, which
scales with the HMF, $\mathrm{d}n/\mathrm{d}\log M_{\mathrm{h}}$, calculated
using the hmf code (Murray et al., 2013) for the fitting function from
Behroozi et al. (2013). A common weight $w=N_{\mathrm{E}}/N_{\mathrm{S}}$ is
assigned to all the halos in the same bin, such that a summation of halo
weights in a given mass bin yields the expected number of halos in the
universe. These weights are then applied to sample galaxies binned in
$M_{\mathrm{UV}}$ to obtain the UVLF, which is essentially a convolution
between the HMF and $M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation including the
full, $M_{\mathrm{h}}$-dependent distribution (see Section 2.4 of Ma et al.,
2018b). Finally, we stress that, compared with Ma et al. (2018b) where only a
subset of the High-Redshift suite was analyzed, we substantially increase the
number of samples of massive halos/bright galaxies in this work (a factor 8
increase of halos with $M_{\mathrm{h}}>10^{10}\,M_{\odot}$ at $z=10$) by
considering the full High-Redshift suite as in Ma et al. (2019), thereby
extending the magnitude down to which the UVLF at $z>10$ can be reliably
determined to $M_{\mathrm{UV}}<-20$, overlapping with the bright-end UVLF
probed by JWST.

##  3 Results

Figure 1: Top: UV magnitude–halo mass relations at $z=8$–12. Data for
individual galaxies are denoted by the grey dots (no smoothing applied to the
SFH). The thick solid curves indicate the range of the 5th and 95th
percentiles in the “bursty” and “smoothed” cases, from which the suppression
of bright galaxy number counts due to smoothing is apparent. Bottom: UVLFs at
$z=8$–12 derived from the convolution between the UV magnitude–halo mass
relation and the HMF. Dust-free predictions are shown as solid for both
“bursty” and “smoothed” cases, whereas the dust-attenuated scenario is shown
as dashed for only the “bursty” case (Section 2.2.2) for visual clarity.
Constraints from observations are shown by the data points in black for the
spectroscopically-confirmed-only samples (Harikane et al., 2023a) and in grey
for data sets involving photometric candidates (Oesch et al., 2018; Bowler et
al., 2020; Rojas-Ruiz et al., 2020; Bouwens et al., 2021, 2023; Finkelstein et
al., 2022; Leethochawalit et al., 2022; Castellano et al., 2023; Donnan et
al., 2023; Harikane et al., 2023a; Pérez-González et al., 2023). Cases with
larger and smaller smoothing timescale $\tau_{\mathrm{SF}}$ values than the
fiducial one (100 Myr) are shown at $z=10$ to illustrate the impact of SFH
smoothing on the UVLF.

###  3.1 The $M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ Relation and the UVLF

Following the methods outlined in Sections 2.2 and 2.3, we first use our
samples of simulated galaxies to quantify the
$M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation in different redshift regimes,
assuming either “bursty” or “smoothed” SFH. A comparison of the
$M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relations at $z=8$, 10, and 12 from our
simulations is shown in the top row of Figure 1. Overall, galaxies become more
UV-bright at higher $M_{\mathrm{h}}$ and, at a given $M_{\mathrm{h}}$,
$M_{\mathrm{UV}}$ decreases modestly with increasing redshift as a result of
more rapid halo growth at higher redshift. A significant scatter in
$M_{\mathrm{UV}}$ around the median relation that gradually increases towards
lower masses exists, which is a sign of increasing star formation burstiness
at low masses, given the proportionality between $L_{\mathrm{UV}}$ and the
SFR. At a fixed $M_{\mathrm{h}}$, we find a modest trend for the scatter in
$M_{\mathrm{UV}}$ to decrease with decreasing redshift that continues to $z<8$
(not shown). This tentative evidence for the redshift evolution of the UV
variability might be testable using comparisons of SFR indicators sensitive to
different star formation timescales or high-precision measurements of the
halo–galaxy connection with galaxy clustering (see Section 4). From the
comparison between the “bursty” and “smoothed” cases shown by the 5–95th
percentiles (especially in the top middle panel where three “smoothed” cases
with varying $\tau_{\mathrm{SF}}$ are shown), it can be seen that evaluating
$M_{\mathrm{UV}}$ from a smoothed SFH leads to a shallower
$M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation with a reduced scatter in
$M_{\mathrm{UV}}$ at higher masses, which effectively suppresses the
population of UV-bright galaxies at a given $M_{\mathrm{h}}$.

In the bottom row of Figure 1, we show the UVLF at $z=8$–12 implied by the
$M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation. From the comparisons against
recent observational constraints and between the two SFH cases, several key
results are immediately apparent. First, in the fiducial, “bursty” SFH
scenario, the predicted UVLFs agree remarkably well with the observational
constraints available. In particular, our $z\gtrsim 10$ predictions lie safely
above the firm lower bounds set by the dust-uncorrected, spectroscopically-
confirmed samples recently compiled by Harikane et al. (2023a), and they are
also broadly consistent with the variety of measurements based on
photometrically selected candidates (see the caption for details)666We have
verified by bootstrapping 1000 times the simulated galaxy samples that the
statistical uncertainty on the UVLF, especially at the bright end, is small
enough that it does not affect the bright-end comparisons of interest to this
study. In the brightest bin, the 1$\sigma$ statistical uncertainties in
$\log\phi$ estimated from bootstrapping are approximately 0.15 dex, 0.15 dex,
and 0.3 dex at $z=8$, 10, and 12, respectively.. Unlike some other theoretical
predictions (e.g., Mason et al., 2023; Yung et al., 2023), for which a clear
tension with the spec-$z$ lower bounds exists without modifications, our
bursty-case predictions do not require any additional tuning of UV variability
or production efficiency to match observations. Despite uncertainties
associated with the treatment of dust, this good agreement implies that the
UVLFs observed by JWST at $z\gtrsim 10$ are consistent with generally “normal”
SFE and UV production efficiency as predicted by the FIRE-2 simulations. As
demonstrated in Ma et al. (2018b), the relation between $M_{*}$ and
$M_{\mathrm{h}}$ in these simulations is broadly consistent with
extrapolations from lower $z$ where empirical analyses show that the SFE is
strongly suppressed by stellar feedback in low-mass halos (Behroozi et al.,
2013; Tacchella et al., 2018).

Table 1: Dust-free UVLFs at $z=8$, 10, and 12 from the simulated galaxies.

$M_{\mathrm{UV}}$ | $\log\phi$ | $M_{\mathrm{UV}}$ | $\log\phi$ | $M_{\mathrm{UV}}$ | $\log\phi$
---|---|---|---|---|---
$z=8$ | $z=10$ | $z=12$
$-10.5$ | $-0.085$ | $-10.25$ | $-0.207$ | $-9.75$ | $-0.234$
$-12.5$ | $-0.570$ | $-12.25$ | $-0.677$ | $-11.75$ | $-0.971$
$-14.5$ | $-1.206$ | $-14.25$ | $-1.242$ | $-13.75$ | $-1.576$
$-16.5$ | $-1.926$ | $-16.25$ | $-2.124$ | $-15.75$ | $-2.200$
$-18.5$ | $-2.815$ | $-18.25$ | $-3.072$ | $-17.75$ | $-3.282$
$-20.5$ | $-3.872$ | $-20.25$ | $-4.344$ | $-19.75$ | $-4.500$
$-22.5$ | $-5.158$ | $-22.25$ | $-5.902$ |  |

  * •

Notes.
$\phi$ values are quoted in units of $\mathrm{mag^{-1}\,Mpc^{-3}}$. See
Equation (1) for analytic fits to the UVLF over $8<z<12$. For reference, in
the two brightest bins, $\phi$ is extracted from a sample of (39, 17), (39,
13), (93, 17) galaxies at $z=8$, 10, and 12, respectively.

Second, in the contrasting, “smoothed” SFH scenario, a clear deficit of UV-
bright galaxies is seen as a result of suppressed up-scattering in
$M_{\mathrm{UV}}$ of low-mass halos when the SFR is averaged over a long
timescale $\tau_{\mathrm{SF}}$. The underestimated abundance of UV-bright
galaxies reveals the important role played by the burstiness of star formation
in determining the number statistics of galaxies at cosmic dawn. As also shown
by the comparison of different $\tau_{\mathrm{SF}}$ values at $z=10$, smoothed
SFHs with $\tau_{\mathrm{SF}}\gtrsim 100\,$Myr result in bright-end UVLFs that
are too steep compared with observations, especially the photometrically
selected samples, for which the bright-end UVLF can be underpredicted at
$>2\sigma$ level in some cases (Castellano et al., 2023; Donnan et al., 2023).
At $z=12$, predictions of the smoothed SFH are in tension with even the most
conservative lower limits derived from only the spectroscopically confirmed
samples (Harikane et al., 2023a). It is therefore clear that the UVLF serves
as a useful probe of the burstiness in the SFH, as been noted in e.g.,
Furlanetto & Mirocha (2022) and Shen et al. (2023), although in practice it
can be challenging to extract the burstiness information from only the UVLF
measurements (see Section 4). The overall shallower
$M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation when the SFH is smoothed also
leads to slightly steeper slope at the faint end, although the effect is much
smaller than the suppression at the bright end.

Figure 2: Dust-free UVLFs at $z=8$, 10, and 12 predicted by the FIRE-2
simulations and from the literature. The binned and the best-fit, double-power
law UVLFs are denoted by the crosses and solid curves, as specified in Table 1
and Equation (1), respectively. Several example dust-free predictions from
other cosmological hydrodynamical simulations, including MillenniumTNG
(dashed, Kannan et al., 2023), FLARES (dotted, Vijayan et al., 2021; Wilkins
et al., 2023), and CoDa II (dotted and only at $z=8$ and 10, Ocvirk et al.,
2020; Dawoodbhoy et al., 2023) are also plotted for comparison.

The binned UVLFs without dust attenuation extracted from our simulations at
$z=8$, 10, and 12 are summarized in Table 1. As has been demonstrated in
Figure 1, dust attenuation only modestly affects the UVLF at the very bright
end, reducing $\phi$ (in the brightest bin) by approximately 0.4, 0.25, and
0.01 dex at $z=8$, 10, and 12, respectively. The binning scheme is chosen such
that the brightest $M_{\mathrm{UV}}$ bin contains more than ten simulated
galaxies for robust statistics. Meanwhile, we fit the dust-free UVLF at $8\leq
z\leq 12$ assuming a universal double-power law (DPL) in $M_{\mathrm{UV}}$,

$\Phi(M_{\mathrm{UV}})=\frac{0.4(\ln
10)\,10^{\phi_{*}}}{10^{0.4(\alpha+1)(M^{*}_{\mathrm{UV}}-M_{\mathrm{UV}})}+10^{0.4(\beta+1)(M^{*}_{\mathrm{UV}}-M_{\mathrm{UV}})}}.$
(1)

We specify the redshift-dependent DPL parameters $\phi_{*}$,
$M^{*}_{\mathrm{UV}}$, $\alpha$, and $\beta$ in the form of a single power law
as $\phi_{*}(z)=\phi_{*,0}[(1+z)/10]^{\phi_{*,1}}$,
$M^{*}_{\mathrm{UV}}(z)=M^{*,0}_{\mathrm{UV}}[(1+z)/10]^{M^{*,1}_{\mathrm{UV}}}$,
$\alpha_{*}(z)=\alpha_{*,0}[(1+z)/10]^{\alpha_{*,1}}$, and
$\beta_{*}(z)=\beta_{*,0}[(1+z)/10]^{\beta_{*,1}}$, where the best-fit
parameters are found to be $\phi_{*,0}=-2.01$, $\phi_{*,1}=0.68$,
$M^{*,0}_{\mathrm{UV}}=-17.26$, $M^{*,1}_{\mathrm{UV}}=-0.08$,
$\alpha_{*,0}=-0.31$, $\alpha_{*,1}=-0.93$, $\beta_{*,0}=0.68$, and
$\alpha_{*,1}=0.93$.

Figure 2 shows a comparison between the binned and best-fit UVLFs predicted by
our simulations and other theoretical predictions in the literature based on
cosmological hydrodynamical simulations (Ocvirk et al., 2020; Vijayan et al.,
2021; Dawoodbhoy et al., 2023; Kannan et al., 2023; Wilkins et al., 2023).
Overall, our predicted UVLFs show a weaker redshift evolution beyond $z=8$
compared with the predictions from the MillenniumTNG (Kannan et al., 2023) and
CoDa II (Ocvirk et al., 2020; Dawoodbhoy et al., 2023) simulations, which
results in a higher abundance of bright ($M_{\mathrm{UV}}\lesssim-20$)
galaxies at $z\gtrsim 10$. Our bright-end predictions are generally comparable
to those from the FLARES simulations (Vijayan et al., 2021; Wilkins et al.,
2023) in both normalization and slope, despite the vastly different nature of
the simulations and methods to evaluate the UVLF. It is noteworthy, though,
that the FIRE-2 simulations analyzed in this work have significantly higher
resolution ($m_{\mathrm{b}}\approx 7\times 10^{3}\,M_{\odot}$ in FIRE-2 vs.
$m_{\mathrm{b}}\approx 2\times 10^{6}\,M_{\odot}$ in FLARES), which allows us
to predict the UVLFs at $8\leq z\leq 12$ down to $M_{\mathrm{UV}}\sim-10$ vs.
the FLARES predictions down to $M_{\mathrm{UV}}\sim-18$. We have also verified
that UVLFs in this work and from Ma et al. (2018b, 2019) are in good agreement
in the overlapping regime.

###  3.2 UV Luminosity Density

By integrating the predicted UVLFs, we can derive the UV luminosity density,
$\rho_{\mathrm{UV}}$, as a function of time, which traces the cosmic star
formation rate density (SFRD). Since at $z\gtrsim 10$ only the brightest end
($M_{\mathrm{UV}}\ll M_{\mathrm{UV,*}}$) of the UVLF has been probed, we
follow Harikane et al. (2023a) to compare the UV luminosity density
contributed by galaxies brighter than $M_{\mathrm{UV}}=-18$, namely
$\rho_{\mathrm{UV,bright}}=\rho_{\mathrm{UV}}(M_{\mathrm{UV}}<-18)$, which
corresponds to the contribution from halos with $M_{\mathrm{h}}\gtrsim
10^{10}\,M_{\odot}$ at $z=10$. The unconstrained contribution by fainter,
lower-mass galaxies is highly sensitive to the faint-end slope of the UVLF and
might even outweigh $\rho_{\mathrm{UV,bright}}$ (Sun & Furlanetto, 2016), but
the comparison restricted to $M_{\mathrm{UV}}<-18$ galaxies still serves as a
useful test of the overall abundance of bright, massive galaxies and their SFE
at cosmic dawn777Results from this work, Harikane et al. (2018, 2023b), and
Bouwens et al. (2023) are integrated down to $M_{\mathrm{UV,lim}}=-18$,
whereas the rest are down to $M_{\mathrm{UV,lim}}=-17$. Figure 3 thus shows
conservatively that our simulations without smoothing predict enough total UV
emission compared with observations, regardless of the modest difference in
$M_{\mathrm{UV,lim}}$..

Figure 3 shows a comparison of the cumulative UV luminosity density between
the dust-attenuated predictions from our simulations and a compilation of
constraints from observations and theoretical forecasts in the literature.
Throughout, dust-attenuated predictions from models/simulations (curves) are
compared with observations (data points), which are dust-uncorrected. Over
$8\leq z\leq 12$, dust-attenuated luminosity densities predicted by our
simulations without smoothing the SFH are fully consistent with observations
of both photometric galaxy candidates and spectroscopically-confirmed galaxies
that provide firm lower limits. Due to the integrated nature of
$\rho_{\mathrm{UV}}$, the “smoothed” case appears more consistent with the
spec-$z$-only lower limits here than at the bright end of the UVLF as shown in
Figure 2. In both cases with dust attenuation, a power-law evolution of
$\rho_{\mathrm{UV}}\propto(1+z)^{-0.3}$ over $8\leq z\leq 12$ is implied,
which appears more gradual compared with the predictions by some previously
proposed semi-analytic/semi-empirical models, such as in Mason et al. (2015)
and Harikane et al. (2018).

Figure 3: The cumulative UV luminosity density
$\rho_{\mathrm{UV}}(<M_{\mathrm{UV,lim}})$ integrated down to
$M_{\mathrm{UV,lim}}\simeq-18$ with dust attenuation included (see Section
2.2.2). At $z\gtrsim 10$, some theoretical models (e.g., Mason et al., 2015;
Harikane et al., 2018) underestimate $\rho_{\mathrm{UV}}$ compared with
observational constraints based on photometric galaxy candidates (e.g.,
Bouwens et al., 2023; Donnan et al., 2023; McLeod et al., 2023; Pérez-González
et al., 2023) and/or spectroscopically-confirmed galaxies as firm lower limits
(Harikane et al., 2023a). Predictions from our “bursty” case are broadly
consistent with both photometric and spectroscopic samples and show a slightly
weaker redshift evolution $\rho_{\mathrm{UV}}\propto(1+z)^{-0.3}$ over $8\leq
z\leq 12$.

##  4 Discussion and Conclusions

We have demonstrated that the FIRE-2 simulations with a multi-channel
implementation of standard stellar feedback processes can reproduce well the
observed abundance of UV-bright galaxies at $z\gtrsim 10$, including both the
photometrically selected candidates and the spectroscopically confirmed
sources recently discovered by JWST. We further showed that the bursty SFH
predicted to be common in galaxies at cosmic dawn is important for explaining
the bright-end of the UVLF. With burstiness included, the simulations
demonstrate that a boosted UV emissivity due to, e.g., an enhanced SFE, a top-
heavy IMF, AGN contributions, or Population III stars (see e.g., Harikane et
al., 2023b, c), is not necessary to explain the bright-end UVLF at $z\gtrsim
10$. (This is of course not to say that none of these other effects could be
present in the real universe, so it certainly remains interesting to
investigate these other possibilities!) Compared to semi-analytic/empirical
models Mason et al. (2023); Mirocha & Furlanetto (2023); Shen et al. (2023);
Yung et al. (2023), our predictions based on the FIRE-2 simulations avoid ad
hoc fine-tuning of the $M_{\mathrm{UV}}$–$M_{\mathrm{h}}$ relation to match
observations.

Though not shown explicitly in this Letter, we have verified that the stellar
mass–halo mass (SMHM) relation, as a measure of the time-integrated, galaxy-
scale SFE, $f_{\star}\equiv M_{\star}/(f_{\mathrm{b}}M_{\mathrm{h}}$) (where
$f_{\mathrm{b}}=\Omega_{\rm b}/\Omega_{\rm m}$ is the cosmic baryon fraction),
barely evolves over $5\leq z\leq 12$ in our simulations. This is consistent
with the previous results presented in Ma et al. (2018b) based on a subset of
the full High-Redshift suite considered in this work (see their Figure 4).
These results indicate that $f_{\star}$ changes from approximately $10^{-3.3}$
to $10^{-1.5}$ as $M_{\mathrm{h}}$ increases from $10^{8}\,M_{\odot}$ to
$10^{11}\,M_{\odot}$ following a simple power law of slope $\sim 0.6$ in log-
log space. Thus, even though star formation is bursty, the galaxy-scale SFE is
not strongly enhanced in these simulations relative to, e.g., an extrapolation
of the SMHM relation empirically determined at lower redshift (Behroozi et
al., 2019). In particular, our simulations do not appear to realize the
“feedback-free starburst” scenario predicted by Dekel et al. (2023) using
analytic arguments, which would result in $f_{\star}$ values up to order-
unity.888While the FIRE-2 simulations assume a local, instantaneous SFE of
100% per free-fall time, this only applies in dense, self-gravitating gas (see
the methods in Hopkins et al., 2018). On galaxy and molecular cloud scales,
stellar feedback generally regulates the SFE to much lower values (e.g.,
Grudić et al., 2018; Orr et al., 2018; Gurvich et al., 2020).

We note that Pallottini & Ferrara (2023) also recently used a set of
cosmological zoom-in simulations (SERRA; Pallottini et al. 2022) to
investigate some implications of stochastic star formation in early galaxies
for the abundance of $z\gtrsim 10$ galaxies observed by JWST. By
characterizing the distribution of time-dependent variations in the SFR of
individual galaxies, they concluded that the predicted SFR variability cannot
account for the required boost suggested by some recent literature to match
the observed UVLF at $z\gtrsim 10$ (Mirocha & Furlanetto, 2023; Shen et al.,
2023). However, Pallottini & Ferrara (2023) did not self-consistently derive
the UVLF from their simulations. Since other physical factors such as the SFE
also impact the UVLF, in addition to burstiness (Mirocha & Furlanetto, 2023;
Muñoz et al., 2023), in order to unambiguously gauge the importance of bursty
star formation it is desirable to perform a self-consistent, end-to-end study
of the UVLF as we do in this work.

Looking ahead, a detailed characterization of the SFR variability on different
timescales will shed light on the physical processes at play in the build-up
of galaxies at early times, as has been demonstrated in recent work using
periodogram (Pallottini & Ferrara, 2023) or more generally power spectral
density (PSD) analysis (Iyer et al., 2020; Tacchella et al., 2020). Moreover,
various implications of bursty star formation should be explicitly considered
when interpreting observations of high$-z$ galaxies. For example, Sun et al.
(2023a) showed that SFR variability introduces important selection effects in
rest UV-selected samples. Since most galaxies at cosmic dawn may form stars in
a highly bursty manner, the impact of burstiness on galaxy number statistics
also raises questions about how to reliably constrain cosmology with high-$z$
galaxy observations (Sabti et al., 2023). At the same time, it is of great
interest to investigate how to observationally characterize the time
variability of star formation and its mass and redshift dependence, e.g. using
SFR indicators sensitive to different timescales (Sparre et al., 2017; Flores
Velázquez et al., 2021; Sun et al., 2023b) or the spatial clustering of
galaxies (Muñoz et al., 2023). Quantifying the effects of bursty star
formation on statistics such as galaxy clustering is a critical stepping stone
towards the usage of high-$z$ galaxies as robust cosmological probes.

The authors thank the anonymous reviewer for comments that helped improve this
Letter, as well as Pratik Gandhi, Yuichi Harikane, and Julian Muñoz for
helpful discussion. GS was supported by a CIERA Postdoctoral Fellowship. CAFG
was supported by NSF through grants AST-2108230 and CAREER award AST-1652522;
by NASA through grants 17-ATP17-0067 and 21-ATP21-0036; by STScI through grant
HST-GO-16730.016-A; and by CXO through grant TM2-23005X. The Flatiron
Institute is supported by the Simons Foundation. AW received support from: NSF
via CAREER award AST-2045928 and grant AST-2107772; NASA ATP grant
80NSSC20K0513; and HST grants AR-15809, GO-15902, GO-16273 from STScI. The
simulations used in this Letter were run on XSEDE computational resources
(allocations TG-AST120025, TG-AST130039, TG-AST140023, and TG-AST140064).
Additional analysis was done using the Quest computing cluster at Northwestern
University.

##  Appendix A Forming the halo/galaxy sample

Throughout, we analyze snapshots of a galaxy in a $\Delta z=0.5$ bin multiple
times per the cadence at which snapshots are stored (every 10–20 Myr). While
the same galaxy from neighbouring snapshots are not strictly independent as
far as $M_{\mathrm{UV}}$ is considered, this method is useful because the
highly time-variable SFR limits the temporal correlation between consecutive
snapshots. It yields a large statistical sample appropriate for UVLF analysis
(see Figure 4) and the sampling cadence does not bias the results, as have
been shown by analyses that randomly exclude approximately half of the samples
(Ma et al., 2018b). At $z=8$, 10, and 12, the UVLF is evaluated from a sample
of approximately 12,000, 9,000, and 4,000 galaxies, respectively. Summing over
the three redshift bins, this yields $\approx 25,000$ galaxy samples in total.

Figure 4: The number of galaxies sampled in each halo mass bin for the 26
simulations inspected for evaluating the UVLF at $z=10\pm 0.5$. This amounts
to a total of $\approx 9000$ galaxies sampled from snapshots of the 26 zoom-in
regions over $9.5<z<10.5$. Simulation IDs are listed in the legend (c.f., Sun
et al., 2023a).

## References

  * Akins et al. (2023) Akins, H. B., Casey, C. M., Allen, N., et al. 2023, arXiv e-prints, arXiv:2304.12347, doi: 10.48550/arXiv.2304.12347
  * Arrabal Haro et al. (2023) Arrabal Haro, P., Dickinson, M., Finkelstein, S. L., et al. 2023, arXiv e-prints, arXiv:2303.15431, doi: 10.48550/arXiv.2303.15431
  * Bakx et al. (2023) Bakx, T. J. L. C., Zavala, J. A., Mitsuhashi, I., et al. 2023, MNRAS, 519, 5076, doi: 10.1093/mnras/stac3723
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57, doi: 10.1088/0004-637X/770/1/57
  * Biagetti et al. (2023) Biagetti, M., Franciolini, G., & Riotto, A. 2023, ApJ, 944, 113, doi: 10.3847/1538-4357/acb5ea
  * Bird et al. (2023) Bird, S., Chang, C.-F., Cui, Y., & Yang, D. 2023, arXiv e-prints, arXiv:2307.10302.  https://arxiv.org/abs/2307.10302
  * Bouwens et al. (2021) Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, AJ, 162, 47, doi: 10.3847/1538-3881/abf83e
  * Bouwens et al. (2023) Bouwens, R. J., Stefanon, M., Brammer, G., et al. 2023, MNRAS, 523, 1036, doi: 10.1093/mnras/stad1145
  * Bowler et al. (2020) Bowler, R. A. A., Jarvis, M. J., Dunlop, J. S., et al. 2020, MNRAS, 493, 2059, doi: 10.1093/mnras/staa313
  * Boylan-Kolchin (2023) Boylan-Kolchin, M. 2023, Nature Astronomy, doi: 10.1038/s41550-023-01937-7
  * Byler et al. (2017) Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44, doi: 10.3847/1538-4357/aa6c66
  * Byrne et al. (2023) Byrne, L., Faucher-Giguère, C.-A., Stern, J., et al. 2023, MNRAS, 520, 722, doi: 10.1093/mnras/stad171
  * Castellano et al. (2023) Castellano, M., Fontana, A., Treu, T., et al. 2023, ApJ, 948, L14, doi: 10.3847/2041-8213/accea5
  * Cochrane et al. (2022) Cochrane, R. K., Hayward, C. C., & Anglés-Alcázar, D. 2022, ApJ, 939, L27, doi: 10.3847/2041-8213/ac951d
  * Cochrane et al. (2019) Cochrane, R. K., Hayward, C. C., Anglés-Alcázar, D., et al. 2019, MNRAS, 488, 1779, doi: 10.1093/mnras/stz1736
  * Cullen et al. (2023) Cullen, F., McLure, R. J., McLeod, D. J., et al. 2023, MNRAS, 520, 14, doi: 10.1093/mnras/stad073
  * Curtis-Lake et al. (2023) Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, Nature Astronomy, 7, 622, doi: 10.1038/s41550-023-01918-w
  * Dawoodbhoy et al. (2023) Dawoodbhoy, T., Shapiro, P. R., Ocvirk, P., et al. 2023, arXiv e-prints, arXiv:2302.08523, doi: 10.48550/arXiv.2302.08523
  * Dayal & Giri (2023) Dayal, P., & Giri, S. K. 2023, arXiv e-prints, arXiv:2303.14239, doi: 10.48550/arXiv.2303.14239
  * Dekel et al. (2023) Dekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., & Li, Z. 2023, MNRAS, doi: 10.1093/mnras/stad1557
  * Domínguez et al. (2015) Domínguez, A., Siana, B., Brooks, A. M., et al. 2015, MNRAS, 451, 839, doi: 10.1093/mnras/stv1001
  * Donnan et al. (2023) Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011, doi: 10.1093/mnras/stac3472
  * Dressler et al. (2023) Dressler, A., Rieke, M., Eisenstein, D., et al. 2023, arXiv e-prints, arXiv:2306.02469, doi: 10.48550/arXiv.2306.02469
  * Eldridge et al. (2017) Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058, doi: 10.1017/pasa.2017.51
  * Endsley et al. (2023) Endsley, R., Stark, D. P., Whitler, L., et al. 2023, arXiv e-prints, arXiv:2306.05295, doi: 10.48550/arXiv.2306.05295
  * Faucher-Giguère (2020) Faucher-Giguère, C.-A. 2020, MNRAS, 493, 1614, doi: 10.1093/mnras/staa302
  * Faucher-Giguère et al. (2009) Faucher-Giguère, C.-A., Lidz, A., Zaldarriaga, M., & Hernquist, L. 2009, ApJ, 703, 1416, doi: 10.1088/0004-637X/703/2/1416
  * Ferrara et al. (2023) Ferrara, A., Pallottini, A., & Dayal, P. 2023, MNRAS, 522, 3986, doi: 10.1093/mnras/stad1095
  * Finkelstein et al. (2022) Finkelstein, S. L., Bagley, M. B., Haro, P. A., et al. 2022, ApJ, 940, L55, doi: 10.3847/2041-8213/ac966e
  * Flores Velázquez et al. (2021) Flores Velázquez, J. A., Gurvich, A. B., Faucher-Giguère, C.-A., et al. 2021, MNRAS, 501, 4812, doi: 10.1093/mnras/staa3893
  * Fudamoto et al. (2020) Fudamoto, Y., Oesch, P. A., Faisst, A., et al. 2020, A&A, 643, A4, doi: 10.1051/0004-6361/202038163
  * Fudamoto et al. (2021) Fudamoto, Y., Oesch, P. A., Schouws, S., et al. 2021, Nature, 597, 489, doi: 10.1038/s41586-021-03846-z
  * Furlanetto & Mirocha (2022) Furlanetto, S. R., & Mirocha, J. 2022, MNRAS, 511, 3895, doi: 10.1093/mnras/stac310
  * Furlanetto & Mirocha (2023) —. 2023, MNRAS, 523, 5274, doi: 10.1093/mnras/stad1799
  * Gnedin (2000) Gnedin, N. Y. 2000, ApJ, 542, 535, doi: 10.1086/317042
  * Gong et al. (2023) Gong, Y., Yue, B., Cao, Y., & Chen, X. 2023, ApJ, 947, 28, doi: 10.3847/1538-4357/acc109
  * Grudić et al. (2018) Grudić, M. Y., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2018, MNRAS, 475, 3511, doi: 10.1093/mnras/sty035
  * Gurvich et al. (2020) Gurvich, A. B., Faucher-Giguère, C.-A., Richings, A. J., et al. 2020, MNRAS, 498, 3664, doi: 10.1093/mnras/staa2578
  * Gurvich et al. (2023) Gurvich, A. B., Stern, J., Faucher-Giguère, C.-A., et al. 2023, MNRAS, 519, 2598, doi: 10.1093/mnras/stac3712
  * Harikane et al. (2023a) Harikane, Y., Nakajima, K., Ouchi, M., et al. 2023a, arXiv e-prints, arXiv:2304.06658, doi: 10.48550/arXiv.2304.06658
  * Harikane et al. (2018) Harikane, Y., Ouchi, M., Ono, Y., et al. 2018, PASJ, 70, S11, doi: 10.1093/pasj/psx097
  * Harikane et al. (2023b) Harikane, Y., Ouchi, M., Oguri, M., et al. 2023b, ApJS, 265, 5, doi: 10.3847/1538-4365/acaaa9
  * Harikane et al. (2023c) Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023c, arXiv e-prints, arXiv:2303.11946, doi: 10.48550/arXiv.2303.11946
  * Hassan et al. (2023) Hassan, S., Lovell, C. C., Madau, P., et al. 2023, arXiv e-prints, arXiv:2305.02703, doi: 10.48550/arXiv.2305.02703
  * Hirano & Yoshida (2023) Hirano, S., & Yoshida, N. 2023, arXiv e-prints, arXiv:2306.11993, doi: 10.48550/arXiv.2306.11993
  * Hopkins (2015) Hopkins, P. F. 2015, MNRAS, 450, 53, doi: 10.1093/mnras/stv195
  * Hopkins et al. (2014) Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581, doi: 10.1093/mnras/stu1738
  * Hopkins et al. (2018) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * Hopkins et al. (2023) Hopkins, P. F., Gurvich, A. B., Shen, X., et al. 2023, MNRAS, doi: 10.1093/mnras/stad1902
  * Hu et al. (2023) Hu, C.-Y., Smith, M. C., Teyssier, R., et al. 2023, ApJ, 950, 132, doi: 10.3847/1538-4357/accf9e
  * Inayoshi et al. (2022) Inayoshi, K., Harikane, Y., Inoue, A. K., Li, W., & Ho, L. C. 2022, ApJ, 938, L10, doi: 10.3847/2041-8213/ac9310
  * Iyer et al. (2020) Iyer, K. G., Tacchella, S., Genel, S., et al. 2020, MNRAS, 498, 430, doi: 10.1093/mnras/staa2150
  * Kannan et al. (2023) Kannan, R., Springel, V., Hernquist, L., et al. 2023, MNRAS, 524, 2594, doi: 10.1093/mnras/stac3743
  * Keller et al. (2023) Keller, B. W., Munshi, F., Trebitsch, M., & Tremmel, M. 2023, ApJ, 943, L28, doi: 10.3847/2041-8213/acb148
  * Knollmann & Knebe (2009) Knollmann, S. R., & Knebe, A. 2009, ApJS, 182, 608, doi: 10.1088/0067-0049/182/2/608
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Leethochawalit et al. (2022) Leethochawalit, N., Roberts-Borsani, G., Morishita, T., Trenti, M., & Treu, T. 2022, arXiv e-prints, arXiv:2205.15388, doi: 10.48550/arXiv.2205.15388
  * Looser et al. (2023a) Looser, T. J., D’Eugenio, F., Maiolino, R., et al. 2023a, arXiv e-prints, arXiv:2302.14155, doi: 10.48550/arXiv.2302.14155
  * Looser et al. (2023b) —. 2023b, arXiv e-prints, arXiv:2306.02470, doi: 10.48550/arXiv.2306.02470
  * Ma et al. (2018a) Ma, X., Hopkins, P. F., Boylan-Kolchin, M., et al. 2018a, MNRAS, 477, 219, doi: 10.1093/mnras/sty684
  * Ma et al. (2018b) Ma, X., Hopkins, P. F., Garrison-Kimmel, S., et al. 2018b, MNRAS, 478, 1694, doi: 10.1093/mnras/sty1024
  * Ma et al. (2019) Ma, X., Hayward, C. C., Casey, C. M., et al. 2019, MNRAS, 487, 1844, doi: 10.1093/mnras/stz1324
  * Mason et al. (2015) Mason, C. A., Trenti, M., & Treu, T. 2015, ApJ, 813, 21, doi: 10.1088/0004-637X/813/1/21
  * Mason et al. (2023) —. 2023, MNRAS, 521, 497, doi: 10.1093/mnras/stad035
  * McCaffrey et al. (2023) McCaffrey, J., Hardin, S., Wise, J., & Regan, J. 2023, arXiv e-prints, arXiv:2304.13755, doi: 10.48550/arXiv.2304.13755
  * McLeod et al. (2023) McLeod, D. J., Donnan, C. T., McLure, R. J., et al. 2023, arXiv e-prints, arXiv:2304.14469, doi: 10.48550/arXiv.2304.14469
  * Melia (2023) Melia, F. 2023, MNRAS, 521, L85, doi: 10.1093/mnrasl/slad025
  * Mirocha & Furlanetto (2023) Mirocha, J., & Furlanetto, S. R. 2023, MNRAS, 519, 843, doi: 10.1093/mnras/stac3578
  * Muñoz et al. (2023) Muñoz, J. B., Mirocha, J., Furlanetto, S., & Sabti, N. 2023, arXiv e-prints, arXiv:2306.09403, doi: 10.48550/arXiv.2306.09403
  * Muratov et al. (2015) Muratov, A. L., Kereš, D., Faucher-Giguère, C.-A., et al. 2015, MNRAS, 454, 2691, doi: 10.1093/mnras/stv2126
  * Murray et al. (2013) Murray, S. G., Power, C., & Robotham, A. S. G. 2013, Astronomy and Computing, 3, 23, doi: 10.1016/j.ascom.2013.11.001
  * Naidu et al. (2022a) Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022a, ApJ, 940, L14, doi: 10.3847/2041-8213/ac9b22
  * Naidu et al. (2022b) Naidu, R. P., Oesch, P. A., Setton, D. J., et al. 2022b, arXiv e-prints, arXiv:2208.02794, doi: 10.48550/arXiv.2208.02794
  * Noh & McQuinn (2014) Noh, Y., & McQuinn, M. 2014, MNRAS, 444, 503, doi: 10.1093/mnras/stu1412
  * Ocvirk et al. (2020) Ocvirk, P., Aubert, D., Sorce, J. G., et al. 2020, MNRAS, 496, 4087, doi: 10.1093/mnras/staa1266
  * Oesch et al. (2018) Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105, doi: 10.3847/1538-4357/aab03f
  * Oke & Gunn (1983) Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713, doi: 10.1086/160817
  * Orr et al. (2018) Orr, M. E., Hayward, C. C., Hopkins, P. F., et al. 2018, MNRAS, 478, 3653, doi: 10.1093/mnras/sty1241
  * Padmanabhan & Loeb (2023) Padmanabhan, H., & Loeb, A. 2023, arXiv e-prints, arXiv:2306.04684, doi: 10.48550/arXiv.2306.04684
  * Pallottini & Ferrara (2023) Pallottini, A., & Ferrara, A. 2023, arXiv e-prints, arXiv:2307.03219, doi: 10.48550/arXiv.2307.03219
  * Pallottini et al. (2022) Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621, doi: 10.1093/mnras/stac1281
  * Parashari & Laha (2023) Parashari, P., & Laha, R. 2023, arXiv e-prints, arXiv:2305.00999, doi: 10.48550/arXiv.2305.00999
  * Pérez-González et al. (2023) Pérez-González, P. G., Costantin, L., Langeroodi, D., et al. 2023, arXiv e-prints, arXiv:2302.02429, doi: 10.48550/arXiv.2302.02429
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Qin et al. (2023) Qin, Y., Balu, S., & Wyithe, J. S. B. 2023, arXiv e-prints, arXiv:2305.17959, doi: 10.48550/arXiv.2305.17959
  * Reddy et al. (2018) Reddy, N. A., Oesch, P. A., Bouwens, R. J., et al. 2018, ApJ, 853, 56, doi: 10.3847/1538-4357/aaa3e7
  * Robertson (2022) Robertson, B. E. 2022, ARA&A, 60, 121, doi: 10.1146/annurev-astro-120221-044656
  * Robertson et al. (2023) Robertson, B. E., Tacchella, S., Johnson, B. D., et al. 2023, Nature Astronomy, 7, 611, doi: 10.1038/s41550-023-01921-1
  * Rojas-Ruiz et al. (2020) Rojas-Ruiz, S., Finkelstein, S. L., Bagley, M. B., et al. 2020, ApJ, 891, 146, doi: 10.3847/1538-4357/ab7659
  * Sabti et al. (2023) Sabti, N., Muñoz, J. B., & Kamionkowski, M. 2023, arXiv e-prints, arXiv:2305.07049, doi: 10.48550/arXiv.2305.07049
  * Shen et al. (2023) Shen, X., Vogelsberger, M., Boylan-Kolchin, M., Tacchella, S., & Kannan, R. 2023, arXiv e-prints, arXiv:2305.05679, doi: 10.48550/arXiv.2305.05679
  * Shen et al. (2022) Shen, X., Vogelsberger, M., Nelson, D., et al. 2022, MNRAS, 510, 5560, doi: 10.1093/mnras/stab3794
  * Sipple & Lidz (2023) Sipple, J., & Lidz, A. 2023, arXiv e-prints, arXiv:2306.12087, doi: 10.48550/arXiv.2306.12087
  * Sparre et al. (2017) Sparre, M., Hayward, C. C., Feldmann, R., et al. 2017, MNRAS, 466, 88, doi: 10.1093/mnras/stw3011
  * Stern et al. (2021) Stern, J., Faucher-Giguère, C.-A., Fielding, D., et al. 2021, ApJ, 911, 88, doi: 10.3847/1538-4357/abd776
  * Sun et al. (2023a) Sun, G., Faucher-Giguère, C.-A., Hayward, C. C., & Shen, X. 2023a, arXiv e-prints, arXiv:2305.02713, doi: 10.48550/arXiv.2305.02713
  * Sun & Furlanetto (2016) Sun, G., & Furlanetto, S. R. 2016, MNRAS, 460, 417, doi: 10.1093/mnras/stw980
  * Sun et al. (2023b) Sun, G., Lidz, A., Faisst, A. L., & Faucher-Giguère, C.-A. 2023b, MNRAS, doi: 10.1093/mnras/stad2000
  * Tacchella et al. (2018) Tacchella, S., Bose, S., Conroy, C., Eisenstein, D. J., & Johnson, B. D. 2018, ApJ, 868, 92, doi: 10.3847/1538-4357/aae8e0
  * Tacchella et al. (2020) Tacchella, S., Forbes, J. C., & Caplar, N. 2020, MNRAS, 497, 698, doi: 10.1093/mnras/staa1838
  * Tamura et al. (2019) Tamura, Y., Mawatari, K., Hashimoto, T., et al. 2019, ApJ, 874, 27, doi: 10.3847/1538-4357/ab0374
  * Topping et al. (2023) Topping, M. W., Stark, D. P., Endsley, R., et al. 2023, arXiv e-prints, arXiv:2307.08835, doi: 10.48550/arXiv.2307.08835
  * Trinca et al. (2023) Trinca, A., Schneider, R., Valiante, R., et al. 2023, arXiv e-prints, arXiv:2305.04944, doi: 10.48550/arXiv.2305.04944
  * Vijayan et al. (2021) Vijayan, A. P., Lovell, C. C., Wilkins, S. M., et al. 2021, MNRAS, 501, 3289, doi: 10.1093/mnras/staa3715
  * Vogelsberger et al. (2020) Vogelsberger, M., Nelson, D., Pillepich, A., et al. 2020, MNRAS, 492, 5167, doi: 10.1093/mnras/staa137
  * Wetzel & Garrison-Kimmel (2020) Wetzel, A., & Garrison-Kimmel, S. 2020, GizmoAnalysis: Read and analyze Gizmo simulations, Astrophysics Source Code Library, record ascl:2002.015.  http://ascl.net/2002.015
  * Wetzel et al. (2023) Wetzel, A., Hayward, C. C., Sanderson, R. E., et al. 2023, ApJS, 265, 44, doi: 10.3847/1538-4365/acb99a
  * Wetzel et al. (2016) Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, ApJ, 827, L23, doi: 10.3847/2041-8205/827/2/L23
  * Wilkins et al. (2023) Wilkins, S. M., Vijayan, A. P., Lovell, C. C., et al. 2023, MNRAS, 519, 3118, doi: 10.1093/mnras/stac3280
  * Witstok et al. (2023) Witstok, J., Shivaei, I., Smit, R., et al. 2023, arXiv e-prints, arXiv:2302.05468, doi: 10.48550/arXiv.2302.05468
  * Yan et al. (2023) Yan, H., Ma, Z., Ling, C., Cheng, C., & Huang, J.-S. 2023, ApJ, 942, L9, doi: 10.3847/2041-8213/aca80c
  * Yung et al. (2023) Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Wilkins, S. M., & Gardner, J. P. 2023, arXiv e-prints, arXiv:2304.04348, doi: 10.48550/arXiv.2304.04348
  * Zavala et al. (2023) Zavala, J. A., Buat, V., Casey, C. M., et al. 2023, ApJ, 943, L9, doi: 10.3847/2041-8213/acacfe


# Magnetic field effect on tunneling through triple barrier in AB bilayer
graphene

Mouhamadou Hassane Saley<EMAIL_ADDRESS>Laboratory of Theoretical
Physics, Faculty of Sciences, Chouaïb Doukkali University, PO Box 20, 24000 El
Jadida, Morocco  Ahmed Jellal<EMAIL_ADDRESS>Laboratory of Theoretical
Physics, Faculty of Sciences, Chouaïb Doukkali University, PO Box 20, 24000 El
Jadida, Morocco  Canadian Quantum Research Center, 204-3002 32 Ave Vernon, BC
V1T 2L7, Canada

###### Abstract

We investigate electrons tunneling in AB bilayer graphene through a triple
electrostatic barrier in the presence of a perpendicular magnetic field. By
way of the transfer matrix method and using the continuity conditions at the
different interfaces, transmission probability is determined. Additional
resonances appear for two-band tunneling and at normal incidence, and their
number is proportional to the value of $U_{4}$ in cases where $U_{2}<U_{4}$ is
present. However, when $U_{2}>U_{4}$, anti-Klein tunneling increases with
$U_{2}$. The transmission probability exhibits an interesting oscillatory
behavior when $U_{3}>U_{2}=U_{4}$ and $U_{3}<U_{2}=U_{4}$. For fixed energy
$E=0.39\gamma_{1}$, increasing barrier widths increases the number of
oscillations and decreases Klein tunneling. An interlayer bias creates a gap
for: $U_{2}<U_{3}<U_{4}$ and $U_{3}>U_{2}=U_{4}$. In the four-band tunneling
case, transmission decreases in $T^{+}_{+}$, $T^{-}_{+}$ and $T^{-}_{-}$
channels in comparison with the single barrier case. It does, however,
increase for $T^{+}_{-}$ when compared to the single barrier case.
Transmission is suppressed in the gap region when an interlayer bias is
introduced. Our results are relevant for electron confinements in AB bilayer
graphene and for the development of graphene-based transistors.

###### pacs:

72.80.Vp, 73.21.Ac, 73.23.Ad
Keywords: Bilayer graphene, AB-Stacking, triple barrier, magnetic field,
energy spectrum, transmission channels, Klein tunneling

##  I Introduction

Since its discovery in 2004, graphene [1] has been the focus of multiple
theoretical and experimental investigations. This is due to its remarkable
features on the one hand, and its broad range of potential applications in
both fundamental and technological sciences on the other. In fact, unlike
charge carriers in ordinary materials, those in graphene behave like
relativistic massless particles and possess a linear energy dispersion at the
Dirac points [2, 3, 4, 5]. As a result, graphene exhibits unusual electronic
properties such as an unconventional quantum Hall effect [5, 6, 7], the Klein
paradox [2, 8, 9], and a minimum conductivity [6, 10]. Graphene is a promising
material for applications such as carbon-based transistors [11, 12, 13, 14],
optoelectronic devices [15, 16, 17], and strain sensors [18, 19, 20, 21] due
to its high electronic mobility [6, 22, 23], high thermal conductivity [24,
25], and optical and mechanical qualities [26, 27]. However, since the
graphene conduction and valence bands touch each other, an energy gap must be
opened to confine electrons in graphene and enable transistors to switch off
[28, 29]. One of the techniques to open and control a gap is the application
of an external electric field on AB bilayer graphene [30, 31, 32].

AB bilayer graphene results from the stacking of two layers of graphene in
accordance with Bernal’s method [33]. Similar to monolayer graphene, AB
bilayer graphene has several intriguing characteristics. It does, in fact,
display a peculiar quantum Hall effect [34, 35] that is distinct from the Hall
effect seen in monolayer graphene. In contrast to monolayer graphene, electron
tunneling in bilayer graphene is characterized by anti-Klein tunneling [2, 36,
37], i.e., a perfect reflection. Its energy spectrum shows four parabolic
bands [30, 38]. Two of them touch each other at zero energy, and two others
are separated by an energy equal to the interlayer coupling $\gamma_{1}=0.4$
eV [30, 39]. The ability to open and control a gap makes AB bilayer graphene a
major asset for nanoelectronic applications [39, 40].

Many studies on the transport properties of bilayer graphene for two-band
$(E<\gamma_{1})$ and four-band $(E>\gamma_{1})$ tunneling have been published
recently [41, 42, 43, 44, 45, 46]. Motivated by ref. [43], we intend to study
the tunneling of electrons in bilayer graphene through a triple potential
barrier system in the presence of a perpendicular magnetic field. First, we
assessed the case where the energy is inferior to $\gamma_{1}$. It was
discovered that in the triple barrier case and for $U_{2}<U_{4}$, a resonance
appears at normal incidence for energies less than $U_{4}$, which is not the
case for the single barrier [43]. The number of resonances is found to
increase as the value of $U_{4}$ increases. Transmission is zero in the same
region when $U_{2}>U_{4}$, and anti-Klein tunneling increases with $U_{2}$ but
does not change when $U_{3}$ changes. An interesting number of oscillations
have appeared in the triple barrier cases in the energy region greater than
$U_{4}$, particularly when $U_{3}>U_{2}=U_{4}$ and when $U_{3}<U_{2}=U_{4}$.
These oscillation numbers are higher than those obtained in refs. [42, 43, 44,
45, 46]. When we fix the value of energy at $E=0.39\gamma_{1}$ and increase
the barrier widths, Klein-tunneling decreases and the number of oscillations
increases more, even at non normal incidence. When $U_{2}<U_{3}<U_{4}$ and
$U_{3}>U_{2}=U_{4}$ are used, an interlayer bias opens a gap, whereas it has
no effect when $U_{2}>U_{4}$, $U3<U_{2}<U_{4}$, and $U3<U_{2}=U_{4}$ are used.
Subsequently, we have considered the case of energy superior to $\gamma_{1}$.
The transmission decreased in the channels $T^{+}_{+}$, $T^{-}_{+}$ and
$T^{-}_{-}$ compared to the result in a single barrier[43] while in the
$T^{+}_{-}$ channel it increased. Transmission is zero in the gap region for a
non-zero interlayer bias of $(\delta_{3}=0.3\gamma_{1})$, contrary to the
findings in refs. [46, 47, 48].

This work is organized as follows. In Sec. II, we present the theoretical
model and determine the eigenvectors and eigenvalues of the system. Sec. III
deals with the determination of the transmission probability using the
continuity conditions and the transfer matrix method. In Sec. IV, we
numerically present and discuss the obtained results. Finally, we provide a
summary of the results in Sec. V.

##  II Theory and methods

Figure 1: Schematic representation of five regions, including the triple
barrier and magnetic field.

We consider a system of triple electrostatic barriers through which electrons
in bilayer graphene can tunnel. As depicted in Fig. 1, our system is made up
of five regions, each denoted by $j$. In addition to the electric field, the
central region $(j=3)$ is subjected to a magnetic field. The following
Hamiltonian can describe our system [4, 42]:

$H=\begin{pmatrix}V^{+}&v_{F}\pi^{\dagger}&-v_{4}\pi^{\dagger}&v_{3}\pi\\\
v_{F}\pi&V^{+}&\gamma_{1}&-v_{4}\pi^{\dagger}\\\
-v_{4}\pi&\gamma_{1}&V^{-}&v_{F}\pi^{\dagger}\\\
v_{3}\pi^{\dagger}&-v_{4}\pi&v_{F}\pi&V^{-}\end{pmatrix}$ (1)

where $\pi=p_{x}+i(p_{y}+eA_{y}(x))$ is the in-plan momenta, $v_{F}=10^{6}$
m/s is the Fermi velocity. The potentials $V^{+}$ and $V^{-}$ on the first and
second layers, respectively, are defined by

$V_{j}^{\pm}=\begin{cases}0,&j=1,5\\\ U_{j}+\xi\delta_{j},&j=2,3,4\end{cases}$
(2)

such that $\xi=+1$ and $\xi=-1$ correspond to the first and second layers. The
barrier strength is represented by the parameter $U_{j}$, and the interlayer
bias is represented by the parameter $\delta_{j}$. The magnetic field is
applied perpendicularly to the graphene layers, and it can be written in terms
of Heaviside step function $\Theta$

$B(x,y)=B\Theta\left[\left(b-x\right)\left(c-x\right)\right]$ (3)

where $B$ is a constant. The component of the vector potential $A_{y}(x)$ in
the Landau gauge is defined as

$A_{y}(x)=\frac{\hbar}{el_{B}^{2}}\begin{cases}b,&\text{if}\quad x<b\\\
x,&\text{if}\quad b<x<c\\\ c,&\text{if}\quad x>c\end{cases}$ (4)

where $l_{B}=\sqrt{\hbar/eB}$ is the magnetic length. The parameters $v_{3}$
and $v_{4}$ have been shown to have no effect on the band structure at high
energy or on the transmission at low energy [30, 42, 49]. Then, we neglect
them and write the Hamiltonian (1) as

$H_{j}=\begin{pmatrix}U_{j}+\delta_{j}&v_{F}\pi^{\dagger}&0&0\\\
v_{F}\pi&U_{j}+\delta_{j}&\gamma_{1}&0\\\
0&\gamma_{1}&U_{j}-\delta_{j}&v_{F}\pi^{\dagger}\\\
0&0&v_{F}\pi&U_{j}-\delta_{j}\end{pmatrix}$ (5)

in the basis of four component spinor
$\psi^{j}(x,y)=\left[\psi_{A_{1}}^{j},\psi_{B_{1}}^{j},\psi_{A_{2}}^{j},\psi_{B_{2}}^{j}\right]^{\dagger}$
where the symbol $\dagger$ denotes the transpose row vector. As a consequence
of the conservation of momentum $k_{y}$ along the $y$-direction, the spinor
can be written as

$\psi^{j}(x,y)=e^{ik_{y}y}\left[\phi_{A_{1}}^{j},\phi_{B_{1}}^{j},\phi_{A_{2}}^{j},\phi_{B_{2}}^{j}\right]^{\dagger}$
(6)

In order to determine the eigenvalues and the eigenvectors of each region, we
use $H_{j}\psi_{j}=E_{j}\Psi_{j}$ to obtain four coupled differential
equations

$\displaystyle-i\sqrt{2}\vartheta_{0}a\phi_{B_{1}}^{j}$
$\displaystyle=(\varepsilon_{j}-\delta_{j})\phi_{A_{1}}^{j}$ (7a)
$\displaystyle i\sqrt{2}\vartheta_{0}a^{\dagger}\phi_{A_{1}}^{j}$
$\displaystyle=(\varepsilon_{j}-\delta_{j})\phi_{B_{1}}^{j}-\gamma_{1}\phi_{A_{2}}^{j}$
(7b) $\displaystyle-i\sqrt{2}\vartheta_{0}a\phi_{B_{2}}^{j}$
$\displaystyle=(\varepsilon_{j}+\delta_{j})\phi_{A_{2}}^{j}-\gamma_{1}\phi_{B_{1}}^{j}$
(7c) $\displaystyle i\sqrt{2}\vartheta_{0}a^{\dagger}\phi_{A_{2}}^{j}$
$\displaystyle=(\varepsilon_{j}+\delta_{j})\phi_{B_{2}}^{j}$ (7d)

where we have set $\varepsilon_{j}=E_{j}-U_{j}$, $\vartheta_{0}=\frac{\hbar
v_{F}}{lb}$ is the energy scale, $a$ and $a^{\dagger}$ are respectively the
annihilation and creation operators

$\displaystyle
a=\frac{l_{B}}{\sqrt{2}}\left(\partial_{x}+k_{y}+\frac{eA_{y}(x)}{\hbar}\right)$
(8a) $\displaystyle
a^{\dagger}=\frac{l_{B}}{\sqrt{2}}\left(-\partial_{x}+k_{y}+\frac{eA_{y}(x)}{\hbar}\right)$
(8b)

which satisfy the commutation relation
$\left[a,a^{\dagger}\right]=\mathbb{I}$. We eliminate the unknowns from
(7a-7d) step by step and obtain for $\phi_{B_{1}}$

$\left[2\vartheta_{0}^{2}aa^{\dagger}-(\varepsilon_{j}+\delta_{j})^{2}\right]\left[2\vartheta_{0}^{2}a^{\dagger}a-(\varepsilon_{j}-\delta_{j})^{2}\right]\phi_{B_{1}}^{j}=\gamma_{1}^{2}(\varepsilon_{j}^{2}-\delta_{j}^{2})\phi_{B_{1}}^{j}$
(9)

###  II.1 Eigenvalues and eigenvectors in region 3

In region 3, where the vector potential is $A_{y}(x)=\frac{\hbar
x}{el_{B}^{2}}$, we introduce the variable
$z=\sqrt{2}(\frac{x}{l_{B}}+k_{y}l_{B})$ and solve (9) to obtain
$\phi_{B_{1}}$. Then, we substitute the result into (7a-7d) to determine the
rest of the spinor components. For more simplicity, we write the solution in
matrix form

$\psi^{3}(x,y)=G_{3}M_{3}(x)C_{3}$ (10)

where $G_{3}=\mathbb{I}_{4}$ is $4\times 4$ identity matrix,
$C_{3}=\left(c_{+},c_{-},d_{+},d_{-}\right)^{\dagger}$ is a constant, and the
matrix $M_{3}(x)$ is given by

$M_{3}(x)=\begin{pmatrix}\eta_{-}\lambda_{+}\chi_{+,-1}^{+}&\eta_{-}^{*}\lambda_{+}\chi_{-,-1}^{+}&\eta_{-}\lambda_{-}\chi_{+,-1}^{-}&\eta_{-}^{*}\lambda_{-}\chi_{-,-1}^{-}\\\
\chi_{+,0}^{+}&\chi_{-,0}^{+}&\chi_{+,0}^{-}&\chi_{-,0}^{-}\\\
\zeta^{+}\chi_{+,0}^{+}&\zeta^{+}\chi_{-,0}^{+}&\zeta^{-}\chi_{+,0}^{-}&\zeta^{-}\chi_{-,0}^{-}\\\
\eta_{+}^{*}\zeta^{+}\chi_{+,1}^{+}&\eta_{+}\zeta^{+}\chi_{-,1}^{+}&\eta_{+}^{*}\zeta^{-}\chi_{+,1}^{-}&\eta_{+}\zeta^{-}\chi_{-,1}^{-}\end{pmatrix}$
(11)

where $\chi^{\tau}_{\pm,l}=D\left[\lambda_{\tau}\pm l,\pm z\right]$ are the
parabolic cylindrical function with argument $z$ and $l=-1,0,1$. We have set
the quantities

$\displaystyle\lambda_{\tau}=-\frac{1}{2}+\frac{\varepsilon_{3}^{2}+\delta_{3}^{2}}{2\vartheta_{0}^{2}}+\tau\frac{\sqrt{(\vartheta_{0}^{2}-2\varepsilon_{3}\delta_{3})^{2}+\gamma_{1}^{2}(\varepsilon_{3}^{2}-\delta_{3}^{2})}}{2\vartheta_{0}^{2}}$
(12)
$\displaystyle\eta_{\pm}=\frac{-i\sqrt{2}\vartheta_{0}}{\varepsilon_{3}\pm\delta_{3}}$
(13)
$\displaystyle\zeta^{\pm}=\frac{-2\vartheta_{0}^{2}\lambda_{\pm}+(\varepsilon_{3}-\delta_{3})^{2}}{\gamma_{1}(\varepsilon_{3}-\delta_{3})}$
(14)

We solve (12) to obtain the energy spectrum in this region. It is given by

$\varepsilon^{\tau}_{\pm,3}=\pm\frac{1}{\sqrt{6}}\left[\mu^{\frac{1}{3}}+\nu\mu^{\frac{-1}{3}}+2A\right]^{\frac{1}{2}}\\\
+\tau\frac{1}{\sqrt{6}}\left[-6B\sqrt{6}\left(\mu^{\frac{1}{3}}+\nu\mu^{\frac{-1}{3}}+2A\right)^{\frac{-1}{2}}-\left(\mu^{\frac{1}{3}}+\nu\mu^{\frac{-1}{3}}-4A\right)\right]^{\frac{1}{2}}$
(15)

where we have defined the parameters

$\displaystyle\mu=-A^{3}+27B^{2}+9AC+\sqrt{\left(-A^{3}+27B^{2}+9AC\right)^{2}-\nu^{3}}$
(16) $\displaystyle\nu=(A^{2}+3C)$ (17) $\displaystyle
A=\delta_{3}^{2}+(2n+1)\vartheta_{0}^{2}+\frac{\gamma_{1}^{2}}{2}$ (18)
$\displaystyle B=\vartheta_{0}^{2}\delta_{3}$ (19) $\displaystyle
C=\left((2n+1)\vartheta_{0}^{2}-\delta_{3}^{2}\right)^{2}-\vartheta_{0}^{4}+\gamma_{1}^{2}\delta_{3}^{2}$
(20)

and $n=\lambda_{\tau}$ is an integer number.

###  II.2 Eigenvalues and eigenvectors in regions 1, 2, 4, 5

In regions $j=1,2,4,5$, the vector potential is constant and set to be
$A_{y}(x)=\frac{\hbar}{el_{B}^{2}}d_{j}$ where

$d_{j}=\left\\{\begin{array}[]{lll}b,&\text{ if }&x<b\\\ c,&\text{ if
}&x>c\end{array}\right.$ (21)

By solving (9) for $\phi_{B_{1}}$, and substituting the result into (7a-7d),
we obtain a general solution in the matrix form

$\psi^{j}(x,y)=G_{j}M_{j}(x)C_{j}e^{ik_{y}y}$ (22)

where $G_{j}$ and $M_{j}(x)$ are igiven by

$\displaystyle
G_{j}=\begin{pmatrix}f_{-}^{+}&-f_{+}^{+}&f_{-}^{-}&-f_{+}^{-}\\\ 1&1&1&1\\\
h^{+}&h^{+}&h^{-}&h^{-}\\\
h^{+}g_{+}^{+}&-h^{+}g_{-}^{+}&h^{-}g_{+}^{-}&-h^{-}g_{-}^{-}\end{pmatrix}$
(23) $\displaystyle M_{j}(x)=\begin{pmatrix}e^{ik_{j}^{+}x}&0&0&0\\\
0&e^{-ik_{j}^{+}x}&0&0\\\ 0&0&e^{ik_{j}^{-}x}&0\\\
0&0&0&e^{-ik_{j}^{-}x}\end{pmatrix}$ (24)

with the parameters

$\displaystyle f_{\pm}^{\tau}=\hbar v_{F}\frac{k_{j}^{\tau}\pm
i\left(k_{y}+\frac{d_{j}}{l_{B}^{2}}\right)}{\varepsilon_{j}-\delta_{j}}$ (25)
$\displaystyle
h^{\tau}=\frac{\left(\varepsilon_{j}-\delta_{j}\right)^{2}-\left(\hbar
v_{F}\right)^{2}\left[\left(k_{j}^{\tau}\right)^{2}+\left(k_{y}+\frac{d_{j}}{l_{B}^{2}}\right)^{2}\right]}{\left(\varepsilon_{j}-\delta_{j}\right)\gamma_{1}}$
(26) $\displaystyle g_{\pm}^{\tau}=\hbar v_{F}\frac{k_{j}^{\tau}\pm
i\left(k_{y}+\frac{d_{j}}{l_{B}^{2}}\right)}{\varepsilon_{j}+\delta_{j}}$ (27)

The wave vector $k_{j}$ along the $x$-direction is expressed as

$k_{j}^{\tau}=\sqrt{\frac{\varepsilon_{j}^{2}+\delta_{j}^{2}+\tau\sqrt{4\varepsilon_{j}^{2}\delta_{j}^{2}+\gamma_{1}^{2}\left(\varepsilon_{j}^{2}-\delta_{j}^{2}\right)}}{\left(\hbar
v_{F}\right)^{2}}-\left(k_{y}+\frac{d_{j}}{l_{B}^{2}}\right)^{2}}$ (28)

giving rise to the eigenenergies

$\varepsilon_{\pm,j}^{\tau}=\pm\sqrt{\delta_{j}^{2}+\left(\hbar
v_{F}k\right)^{2}+\frac{\gamma_{1}^{2}}{2}+\tau\sqrt{\left(\hbar
v_{F}k\right)^{2}\left(4\delta_{j}^{2}+\gamma_{1}^{2}\right)+\frac{\gamma_{1}^{4}}{4}}}$
(29)

where $k=\left[\left(k_{j}^{\tau}\right)^{2}+k_{y}^{2}\right]^{\frac{1}{2}}$
is the wave vector. It is worth noting that in the incident and transmission
regions $j=1,5$, solutions are obtained by requiring $U_{j}=\delta_{j}=0$.

##  III Transport properties

We will calculate the transmission probability corresponding to the present
system. We do this by imposing continuity conditions at each interface of the
triple barrier structure. Thereafter, we can use the transfer matrix method to
make a connection between the coefficients of the incident region and those of
the transmitted one. These coefficients are given by

$C_{\text{1}}^{\tau}=\begin{pmatrix}\delta_{\tau,1}\\\ r_{+}^{\tau}\\\
\delta_{\tau,-1}\\\ r_{-}^{\tau}\end{pmatrix},\quad
C_{\text{5}}^{\tau}=\begin{pmatrix}t_{+}^{\tau}\\\ 0\\\ t_{-}^{\tau}\\\
0\end{pmatrix}$ (30)

where $\delta_{\tau,\pm 1}$ is the Kronecker delta symbol. The continuity at
interfaces $x=a,b,c,d$ gives rise to

$\displaystyle G_{1}M_{1}(a)C_{1}=G_{2}M_{2}(a)C_{2}$ (31) $\displaystyle
G_{2}M_{2}(b)C_{2}=G_{3}M_{3}(b)C_{3}$ (32) $\displaystyle
G_{3}M_{3}(c)C_{3}=G_{4}M_{4}(c)C_{4}$ (33) $\displaystyle
G_{4}M_{4}(d)C_{5}=G_{5}M_{5}(d)C_{5}.$ (34)

We can now connect the coefficients $C_{\text{1}}^{\tau}$ to
$C_{\text{5}}^{\tau}$

$C_{\text{1}}^{\tau}=NC_{\text{5}}^{\tau}$ (35)

where the matrix transfer $N$ takes the form

$N=\prod_{j=\text{1}}^{\text{4}}M_{j}^{-1}(x_{j})G_{j}^{-1}G_{j+1}M_{j+1}(x_{j})$
(36)

After manipulation, we can write (35) as

$\begin{pmatrix}t_{+}^{\tau}\\\ r_{+}^{\tau}\\\ t_{-}^{\tau}\\\
r_{-}^{\tau}\end{pmatrix}=\begin{pmatrix}N_{11}&0&N_{13}&0\\\
N_{21}&-1&N_{23}&0\\\ N_{31}&0&N_{33}&0\\\
N_{41}&0&N_{43}&-1\end{pmatrix}^{-1}\begin{pmatrix}\delta_{\tau,1}\\\ 0\\\
\delta_{\tau,-1}\\\ 0\end{pmatrix}$ (37)

where $N_{ij}$ are the matrix elements of $N$ (36). Then, we can easily derive
the transmission coefficients from (37). They are given by

$\displaystyle
t_{+}^{\tau}=\frac{N_{33}\delta_{\tau,1}-N_{13}\delta_{\tau,-1}}{N_{11}N_{33}-N_{13}N_{31}}$
(38) $\displaystyle
t_{-}^{\tau}=\frac{N_{11}\delta_{\tau,-1}-N_{31}\delta_{\tau,1}}{N_{11}N_{33}-N_{13}N_{31}}.$
(39)

At this stage, we have to introduce the density of current to find all
transmission channels, which is

$\textbf{j}=v_{F}{\Psi}^{\dagger}\vec{\alpha}\Psi$ (40)

where $\vec{\alpha}$ is a $4\times 4$ matrix with two Pauli matrices
$\sigma_{x}$ on the diagonal, and the rest is zero. Then, using (40), we can
calculate the incident $\textbf{j}_{\text{inc}}$ and transmitted
$\textbf{j}_{\text{tra}}$ current densities. As a result, we obtain the four
transmission channels

$\displaystyle
T_{\pm}^{\tau}=\frac{\left|\textbf{j}_{\text{tra}}\right|}{\left|\textbf{j}_{\text{inc}}\right|}=\frac{k_{5}^{\pm}}{k_{1}^{\tau}}\left|t_{\pm}^{\tau}\right|^{2}$
(41)

##  IV Results and discussions

In the following, we will compute and discuss our numerical results. To begin
with, we will look at the case of energy less than the interlayer coupling
$\gamma_{1}$ (two-band tunneling). Here, only one mode of propagation is
possible, so we have just one transmission channel. Next, we will consider
energy greater than the interlayer coupling $\gamma_{1}$ (four-band
tunneling). In this case, two propagation modes are available, which result in
four transmission channels.

###  IV.1 Two-band Tunneling

Figure 2: (Color online): Transmission as a function of energy $E$, at normal
incidence ($k_{y}=0$), for barrier widths $b_{1}=b_{2}=b_{3}=25$ nm,
$l_{B}=13.5$ nm, and $\delta_{2}=\delta_{3}=\delta_{4}=0$. (a):
$U_{2}<U_{3}<U_{4}$ (green line) and $U_{2}>U_{3}>U_{4}$ (red line). (b): As
in (a) except that $U_{4}=0.5\gamma_{1}$ (green line) and
$U_{2}=0.5\gamma_{1}$ (red line). (c): $U_{3}<U_{2}<U_{4}$ (green line) and
$U_{2}>U_{4}>U_{3}$ (red line). (d): $U_{3}<U_{2}=U_{4}$ (green line) and
$U_{3}>U_{2}=U_{4}$ (red line). The blue line corresponds to the result
obtained in [43] for a single barrier.

The transmission probability is shown in Fig. 2 as a function of incident
energy $E$ at normal incidence $(k_{y}=0)$ for barrier widths of
$b_{1}=b_{2}=b_{3}=25$ nm. The magnetic length is set to $l_{B}=13.5$ nm, as
in the single barrier case discussed in [43], and the interlayer bias is set
to $\delta_{2}=\delta_{3}=\delta_{4}=0$. We plot the transmission for two
configurations of the triple barrier system in Fig. 2(a), and compare it to
the result found in [43] (blue line). In contrast to the single barrier case
(blue line), a transmission resonance occurs for energy lower than $U_{4}$
when $U_{2}<U_{3}<U_{4}$ (green line). In the same energy range, however,
transmission is zero (anti-Klein tunneling [2, 36]) when $U_{2}>U_{3}>U_{4}$
(red line). Transmission oscillates faster in the triple barrier cases (green
and red lines) than in the single barrier (blue line [43]) and double barrier
for zero magnetic field [46] cases for energies greater than $U_{4}$. Fig.
2(b) depicts the transmission with the same parameters as Fig. 2(a), but with
$U_{4}$ (green line) and $U_{2}$ (red line) increased. The number of
resonances increases as $U_{4}$ increases for configuration
$U_{2}<U_{3}<U_{4}$ (green line). In the opposite configuration, we see that
anti-Klein tunneling increases as $U_{2}$ increases. In Fig. 2(c), we reduce
$U_{3}$ so that $U_{3}<U_{2}<U_{4}$ (green line) and $U_{2}>U_{4}>U_{3}$ (red
line) are obtained. Transmission occurs for energies less than $U_{2}$ for
$U_{3}<U_{2}<U_{4}$ (green line), but not for $U_{3}$ (blue line see [43]).
For $U_{2}>U_{4}>U_{3}$ (red line), the transmission is nearly identical to
that shown in Fig. 2(a). We conclude that when $U_{2}>U_{4}$, the transmission
decreases regardless of whether $U_{3}$ is large or small. The transmission
for $U_{3}<U_{2}=U_{4}$ (green line) and $U_{3}>U_{2}=U_{4}$ (red line) is
shown in Fig. 2(d). In both cases, the number of transmission oscillations
increases significantly when compared to the single barrier’s result [43].
However, when $U_{3}>U_{2}=U_{4}$, it is less than when $U_{3}<U_{2}=U_{4}$.
It is also worth noting that such a large number of oscillations are not
observed in [42, 44, 45, 46].

Figure 3: (Color online): The same parameters as in Fig. 2 but now for
$\delta_{3}=0.1\gamma_{1}$ and $\delta_{2}=\delta_{4}=0$.

To investigate the interlayer bias effect, we plot the transmission as a
function of incident energy $E$ in Fig. 3, with the same parameters as in Fig.
2 and $\delta_{3}=0.1\gamma_{1}$ set. The presence of interlayer bias in Fig.
3(a,b) opens a gap for a triple barrier system when $U_{2}<U_{3}<U_{4}$ (green
line), as it does for a single barrier system (blue line [43]). In contrast,
there is no gap when $U_{2}>U_{3}>U_{4}$ (red line), and the transmission
behaves similarly to the result in Fig. 2 with a minor difference. The
interlayer bias has no meaningful effect on transmission in Fig. 3(c) when
$U_{3}<U_{2}<U_{4}$ or when $U_{2}>U_{4}>U_{3}$, no gap is found. However, in
Fig. 3(d), there is a gap when $U_{3}>U_{2}=U_{4}$, but not when
$U_{3}<U_{2}=U_{4}$. Note that, in the triple barrier case, the transmission
does not behave like the results in [43, 44, 45, 46]. For instance, the
transmission is not zero in the region preceding the gap in [43, 44, 45],
while it is zero in the triple barrier case. Moreover, in the double barrier
[46] case, the presence of bound states in the well results in the presence of
resonances within the gap. As can be seen, there is no resonance in the
induced gap in the triple barrier case. Let us also note that in the case of a
triple barrier with zero magnetic field [44], it is necessary to apply an
interlayer bias in at least two regions to obtain a gap, whereas in this case
$\delta_{3}$ can open a gap on its own.

Figure 4: (Color online): Density plot of transmission as a function of the
barrier width ($b_{1}=b_{2}=b_{3}$) and the transverse wave vector $k_{y}$ for
$l_{B}=18.5$ nm and $E=0.39\gamma_{1}$. (a): Single barrier case, (b):
$U_{2}<U_{3}<U_{4}$, (c): $U_{3}<U_{2}<U_{4}$ and (d): $U_{3}<U_{2}=U_{4}$.
The interlayer bias are set at $\delta_{1}=\delta_{2}=\delta_{3}=0.$

In Fig. 4 we show the transmission as a function of barrier width
($b_{1}=b_{2}=b_{3}$) and the transverse wave vector $k_{y}$. The magnetic
length is set at $l_{B}=18.5$ nm, the interlayer bias at
$\delta_{1}=\delta_{2}=\delta_{3}=0$, and the energy at $E=0.39\gamma_{1}$.
When we compare Fig. 4(a) (single barrier case [43]) to Fig. 4(b)
($U_{2}<U_{3}<U_{4}$), we see that the number of oscillations increases
significantly in the triple barrier case as the barrier widths increase. In
Fig. 4(c) additional transmission resonances appear where $U_{3}<U_{2}<U_{4}$,
resulting in more oscillations than in Fig. 4(b). The number of resonances, on
the other hand, decreases when $U_{3}<U_{2}=U_{4}$, as shown in Fig. 4(d). As
a result, the oscillations are less frequent than in Fig. 4(b). However, it is
important to note that the number of oscillations is greater than that of the
single barrier [43] and double-barrier [46] cases. Klein tunneling is found to
be reduced in the triple barrier cases when compared to the single barrier
[43]. In fact, at $\lvert k_{y}\rvert=0.3$ nm-1, anti-Klein tunneling begins
at $b_{1}=b_{2}=b_{3}=3$ nm in triple barrier cases, and at 8 nm in single
barrier case. In addition, for barrier widths greater than 20 nm, the
transmission narrows more in the triple barrier cases than it does in the
single barrier case. As a consequence, Klein tunneling decreases while anti-
Klein tunneling increases.

###  IV.2 Four-band Tunneling

Figure 5: (Color online): Density plot of transmission as a function of
incident energy $E$ and the transverse wave vector $k_{y}$. (Left): Single
barrier case with $U_{3}=2.5\gamma_{1}$. (Right): Triple barrier case with
$U_{2}=1.5\gamma_{1}$, $U_{3}=2.5\gamma_{1}$ and $U_{4}=3\gamma_{1}$. The
barrier widths are set at $b_{1}=b_{2}=b_{3}=15$ nm, $l_{B}=13.5$ nm and
$\delta_{2}=\delta_{3}=\delta_{4}=0.$

We plot the transmission as a function of incident energy $E$ and the
transverse wave vector $k_{y}$ in Fig. 5 for the single barrier case [43] with
$U_{3}=2.5\gamma_{1}$ (left) and the triple barrier case with
$U_{2}<U_{3}<U_{4}$ (right). The barrier widths are set to
$b_{1}=b_{2}=b_{3}=15$ nm, the magnetic length is set to $l_{B}=13.5$ nm, and
the interlayer bias is set to $\delta_{2}=\delta_{3}=\delta_{4}=0$. The cloak
effect [36, 42] occurs in the $T^{+}_{+}$ channel and at normal incidence in
the energy region $U_{3}-\gamma_{1}<E<U_{3}$ in the single barrier case, and
for a wide range of energy $U_{3}-2\gamma_{1}<E<U_{3}+0.5\gamma_{1}$ in the
triple barrier case. As a result, the transmission resonances decrease for
energies less than $U_{3}$ when compared to the single barrier case. There
are, however, more thin resonances than those obtained in [44, 46, 47].
Furthermore, the transmission exhibits symmetrical behavior with respect to
normal incidence ($k_{y}=0$) in [42, 43, 44, 45, 46], whereas in the triple
barrier it is more pronounced in the region
$U_{3}-2\gamma_{1}<E<U_{3}+0.5\gamma_{1}$ for $k_{y}<0$. Transmission
probability increases in the $T^{+}_{-}$ channel compared to the single
barrier case and approaches unit near $E=U_{3}$ in the case of a triple
barrier (right). In addition, the transmission oscillates, unlike [42, 43,
46]. In the $T^{-}_{+}$ channel, the transmission probability decreases for
the triple barrier case and gets close to zero, as found in [44]. In contrast,
in the case of a single barrier [43], it is different from zero. The
transmission in the $T^{-}_{-}$ channel decreases from $E=U_{3}$ in the single
barrier case to $E=U_{3}-\gamma_{1}$ in the triple barrier case. Therefore,
the number of resonances diminishes in the triple barrier case, but the
transmission probability still remains higher than in [45, 46]. Resonances
appear in the case of a single barrier for energies greater than
$U_{3}+\gamma_{1}$, whereas transmission is zero in the case of a triple
barrier.

Figure 6: (Color online): The same parameters as in Fig. 5, but now for
$\delta_{3}=0.3\gamma_{1}$ and $\delta_{2}=\delta_{4}=0.$

Fig. 6 shows the transmission channels for the same parameters as in Fig. 5,
but with $\delta_{3}=0.3\gamma_{1}$ and $\delta_{2}=\delta_{4}=0$. We observe
that the presence of the interlayer bias results in a gap region in all
transmission channels, as found in the single barrier case [43]. In addition,
no transmission is found in the gap region, in contrast to [46, 47, 48]. The
cloak effect in the $T^{+}_{+}$ channel is larger in the triple barrier case
than in the single barrier case, as shown in Fig. 5. However, in the
$T^{+}_{-}$ channel, it occurs for the same energy region in both the triple
barrier and single barrier cases. Furthermore, we observe two thin resonances
in the $T^{+}_{-}$ channel and triple barrier case that are not present in the
single barrier case [43]. The transmission probability in the $T^{-}_{+}$
channel and triple barrier case becomes non-zero in the region of energies
greater than $U_{4}$, in contrast to the result in Fig. 5, where it is close
to zero. However, it still remains lower than in the case of the single
barrier. In the $T^{-}_{-}$ channel, the transmission in both single barrier
and triple barrier cases is similar to the results obtained in Fig. 5, except
that the number of resonances decreases in the single barrier case.

##  V Conclusion

In this paper, we have investigated the transmission of charge carriers in AB
bilayer graphene through a triple potential barrier and in the presence of a
perpendicular magnetic field. The eigenvalues and eigenvectors for each region
were first determined using the four-band Hamiltonian. Afterwards, we
calculated the transmission probability by applying the continuity conditions
at the different edges of the system and using the transfer matrix method. The
number of resonances increases with the value of $U_{4}$ (when
$U_{2}<U_{3}<U_{4}$ ) in the two-band tunneling case at normal incidence and
for $E<U_{4}$, whereas none are present in the single barrier case. However,
when $U_{2}>U_{3}>U_{4}$, we discovered that anti-Klein tunneling increased as
the value of $U_{2}$ increased. Furthermore, we discovered that anti-Klein
tunneling does not depend on the value of $U_{3}$ when $U_{2}>U_{4}$, but
decreases when $U_{3}<U_{2}<U_{4}$. For energy, $E>U_{4}$, the transmission
probability oscillates rapidly in the triple barrier cases compared to the
single-barrier ones. A large number of oscillations were observed when
$U_{3}<U_{2}=U_{4}$ and when $U_{3}>U_{2}=U_{4}$. We discovered that Klein
tunneling decreases as the barrier widths and number of oscillations increase
at energy $E=0.39\gamma_{1}$. When an interlayer bias is introduced, as in the
single barrier case, a gap region is found in cases $U_{2}<U_{3}<U_{4}$ and
$U_{3}>U2=U_{4}$.

In the four-band tunneling, it was seen that the transmission decreased in the
channels $T^{+}_{+}$, $T^{-}_{+}$ and $T^{-}_{-}$ when compared to the case of
a single barrier. By contrast, in the $T^{+}_{-}$ channel, there is an
increase in the transmission probability in comparison with the single barrier
case. Contrary to [46, 47, 48], we observed total suppression of transmission
in the gap region in the presence of an interlayer bias,
$(\delta_{3}=0.3\gamma_{1})$. To sum up, our findings revealed interesting
oscillatory features that provide good regulation of charge carrier
transmission in AB bilayer graphene. Our results, on the other hand,
demonstrated that in the presence of a magnetic field, the triple barrier
system is capable of confining electrons in AB bilayer graphene as well as
providing a logic state on/off to graphene-based transistors.

## References

  * [1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Science 306, 666 (2004).
  * [2] M. I. Katsnelson, K. S. Novoselov and A. K. Geim, Nat. Phys. 2, 620 (2006).
  * [3] Y. Zhang, Z. Jiang, J. P. Small, M. S. Purewal, Y.-W. Tan, M. Fazlollahi, J. D. Chudow, J. A. Jaszczak, H. L. Stormer and P. Kim, Phys. Rev. Lett. 96, 136806 (2006).
  * [4] A. H. C. Neto, F. Guinea and N. M. R. Peres, Rev. Mod. Phys. 81, 109 (2009).
  * [5] V. P. Gusynin and S. G. Sharapov, Phys. Rev. Lett. 95, 146801 (2005).
  * [6] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos and A. A. Firsov, Nature 438, 197 (2005).
  * [7] Y. Zhang, Y.-W. Tan, H. L. Stormer and P. Kim, Nature 438, 201 (2005).
  * [8] A. F. Young and P. Kim, Nat. Phys. 5, 222 (2009).
  * [9] N. Stander, B. Huard and D. Goldhaber-Gordon, Phys. Rev. Lett. 102, 026807 (2009).
  * [10] J. Tworzydło, B. Trauzettel, M. Titov, A. Rycerz and C. W. J. Beenakker, Phys. Rev. Lett. 96, 246802 (2004).
  * [11] M. C. Lemme, T. J. Echtermeyer, M. Baus and H. Kurz, IEEE Electron Device Lett. 28, 282 (2007).
  * [12] Y.-M. Lin, A. Valdes-Garcia, S.-J. Han, D. B. Farmer, I. Meric, Y. Sun, Y. Wu, C. Dimitrakopoulos, A. Grill, P. Avouris and K. A. Jenkins, Science 332, 1294 (2011).
  * [13] R. Cheng, J. Bai, L. Liao, H. Zhou, Y. Chen, L. Liu, Y.-C. Lin, S. Jiang, Y. Huang and X. Duan, Proc. Nat. Acad. Sci. 109, 11588 (2012).
  * [14] B. Zhan, C. Li , J. Yang, G. Jenkins, W. Huang and X. Dong, Small 10, 4042 (2014).
  * [15] F. Bonaccorso, Z. Sun, T. Hasan and A. C. Ferrari, Nat. Photon 4, 611 (2010).
  * [16] M. Liu, X.Yin, E. Ulin-Avila, B. Geng, T. Zentgraf, L. Ju, F. Wang and X. Zhang, Nature 474, 64 (2011).
  * [17] J. Wang, X. Mu, M. Sun and T. Mu, Appl. Mater. Today 16, 1 (2019).
  * [18] Y. R. Jeong, H. Park, S. W. Jin, S. Y. Hong, S.-S. Lee and J. S. Ha, Adv. Funct. Mater. 25, 4228 (2015).
  * [19] J. Zhao, C. He, R. Yang, Z. Shi, M. Cheng, W. Yang, G. Xie, D. Wang, D. Shi and G. Zhang, Appl. Phys. Lett. 101, 063112 (2012).
  * [20] H. Tian, Y. Shu, Y.-L. Cui, W.-T. Mi, Y. Yang, D. Xie and T.-L. Ren, Nanoscale 6, 699 (2014).
  * [21] S. Chun, Y. Choi and W. Park, Carbon 116, 753 (2017).
  * [22] K. I. Bolotin, K. J. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim and H. L. Stormer, Solid State Commun. 146, 351 (2008).
  * [23] S. V. Morozov, K. S. Novoselov, M. I. Katsnelson, F. Schedin, D. C. Elias, J. A. Jaszczak and A. K. Geim, Phys. Rev. Lett. 100, 016602 (2008).
  * [24] A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao and C. N. Lau, Nano lett. 8, 902 (2008).
  * [25] A. K. Geim, Science 324, 1530 (2009).
  * [26] R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres and A. K. Geim, Science 320, 1308 (2008).
  * [27] C. Lee, X. Wei, J. W. Kysar and J. Hone, Science 321, 385 (2008).
  * [28] T. Tudorovskiy, K. J. A. Reijnders and M. I. Katsnelson, Phys. Scr. T 146, 014010 (2012).
  * [29] K. C. Yung, W. M. Wu, M. P. Pierpoint and F. V. Kusmartsev, Contemp. Phys. 54, 233 (2013).
  * [30] E. McCann, D. S. L. Abergel and V. I. Fal’ko, Solid State Commun. 143, 110 (2007).
  * [31] E. McCann, D. S. L. Abergel and V. I. Fal’ko, Eur. Phys. J.: Spec. Top. 148, 91 (2007).
  * [32] Y. Zhang, T.-T. Tang, C. Girit, Z. Hao, M. C. Martin, A. Zettl, M. F. Crommie, Y. R. Shen and F. Wang, Nature 459, 820 (2009).
  * [33] J. D. Bernal, Proc. R. Soc. A 106, 749 (1924).
  * [34] K. S. Novoselov, E. McCann, S. V. Morozov, V. I. Fal’ko, M. I. Katsnelson, U. Zeitler, D. Jiang, F. Schedin and A. K. Geim, Nat. Phys. 2, 177 (2006).
  * [35] E. McCann and V. I. Fal’ko, Phys. Rev. Lett. 96, 086805 (2006).
  * [36] B. Van Duppen and F. M. Peeters, EPL 102, 27001 (2013).
  * [37] L. Dell’Anna, P Majari and M. R. Setare, J. Phys.: Condens. Matter 30, 415301 (2018).
  * [38] E. McCann and M. Koshino, Rep. Prog. Phys. 76, 056503 (2013).
  * [39] Z. Q. Li, E. A. Henriksen, Z. Jiang, Z. Hao, M. C. Martin, P. Kim, H. L. Stormer and D. N. Basov1, Phys. Rev. Lett. 102, 037403 (2009).
  * [40] J. Sławińska, I. Zasada and Z. Klusek, Phys. Rev. B 81, 155433 (2010).
  * [41] M. Barbier, P. Vasilopoulos, F. M. Peeters and J. M. Pereira Jr, Phys. Rev. B 79, 155402 (2009).
  * [42] B. Van Duppen and F. M. Peeters, Phys. Rev. B 87, 205427 (2013).
  * [43] A. Jellal, I. Redouani and H. Bahlouli, Physica E 72, 149 (2015).
  * [44] M. Hassane Saley, A. El Mouhafid, A. Jellal and A. Siari, Ann. Phys. 534, 2200308 (2022).
  * [45] N. Benlakhouy, A. El Mouhafid and A. Jellal, Physica E 134, 114835 (2021).
  * [46] H. Bahlouli H. M. Abdullah, A. El Mouhafid and A. Jellal, Mater. Res. Express 4, 025009 (2017).
  * [47] W. Lu, W. Li, C. Xu and C. Ye, J. Phys. D: Appl. Phys. 48, 285102 (2015).
  * [48] I. Redouani, A. Jellal and H. Bahlouli, J. Low Temp. Phys. 181, 197 (2015).
  * [49] H. Cheraghchi and F. Adinehvand, J. Phys.: Condens. Matter 26, 015302 (2014).