\section{Introduction}

The spectral energy distributions (SEDs) of  blazars are characterized by two non-thermal peaks \citep{Dermer:1993cz}. The first peak, in the infrared to X-ray energy range, is produced by the synchrotron emission from  relativistic electrons in the jet. The second peak, in between X-ray to $\gamma$-ray energy, is believed to be produced either from the synchrotron Self-Compton (SSC) scattering of high-energy electrons with the low-energy self-produced synchrotron photons in the jet \citep{Maraschi:1992iz,Murase:2011cy,Gao:2012sq} or from the external Compton (EC) scattering  with external sources such as photons from the accretion disk or broad-line regions \citep{Sikora:1994zb,Blazejowski:2000ck}. The model which explains these two peaks is in general referred to as the leptonic model. Depending on the position of the synchrotron peak, the blazars are classified into four categories:  low-energy peaked blazars (LBLs, $\nu_{peak} < 10^{14}\, Hz$), intermediate-energy peaked blazars (IBLs, $\nu_{peak}$ between $10^{14}\, Hz$ and $10^{15}\, Hz$), high energy-peaked blazars (HBLs, $\nu_{peak}$ between $10^{15}\, Hz$ and $10^{17}\, Hz$) \citep{Abdo:2009iq} and extreme energy peaked blazar (EHBL, $ >10^{17}$ Hz) \citep{Costamante:2001pu}. In particular, the  EHBLs are characterized by low luminosity and have limited variabilities which make these objects difficult to detect \citep{Costamante:2017xqg,Acciari:2019jzj}.

So far, only a few EHBLs have been detected in the very high energy (VHE, $>\, 100$ GeV) regime and among them some well known objects are 1ES 0229+200, 1ES 0347-232, RGB J0710+591 and 1ES 1101-232 \citep{Costamante:2017xqg}. Among all the EHBLs observed so far, the highest energy peak frequencies correspond to 1ES 0229+200 with the synchrotron peak at $\nu^{peak}_{syn}\simeq 3.5\times 10^{19}$ Hz and the SSC peak at $\nu^{peak}_{IC}\simeq 1.5\times 10^{27}$ Hz \citep{Kaufmann:2011az} respectively.
The nearest and the extensively studied HBLs Markarian 421 (Mrk 421), Markarian 501 (Mrk 501) and 1ES 1959+650 have also shown EHBL-like behavior with harder TeV spectra in which the synchrotron and SSC peaks are shifted towards higher energy during some flaring episodes \citep{Ahnen:2018mtr,Hayashida:2020wez,Aleksi__2015}. The EHBL-like behavior of Mrk 501 was observed during two episodes of very-high energy flaring events of May - July 2005 and June 2012. Also, during a multiwavelength (MW) campaign of 1ES 1959+650 between 29th April to 21st November 2016, the MAGIC telescopes observed multi-TeV flaring events during the nights of 13th, 14th June and 1st July 2016 when the position of the synchrotron peak was in EHBL range and exhibiting extreme HBL-like behavior.

It is very difficult to explain the shift of the second peak in the one-zone leptonic model. However, it has been shown that by using unrealistic model parameters, such as the large values of minimum electron Lorentz factor and the bulk Lorentz factor \citep{Ahnen:2018mtr,Hayashida:2020wez} as well as a low magnetic field, the shift can be explained. Alternative solutions have also been proposed, such as: a two-zone leptonic model, IC scattering of the electron with the cosmic microwave background, a spine-layer structured jet model and variants of the hadronic model \citep{Boettcher:2008fh,Acciari:2019ntl,Ahnen:2018mtr}.

In previous studies, \cite{Sahu:2019kfd,Sahu:2019scf} have used the photohadronic model to explain the multi-TeV flaring from many HBLs very well.
In all these studies it is shown that the Fermi accelerated high energy protons
interacting with the seed SSC photons in the jet can explain the VHE spectrum very
well. The SSC photon flux $\Phi_{SSC}$ in this region is a perfect power-law given as
$\Phi_{SSC}\propto \epsilon_{\gamma}^{\beta} \propto E_{\gamma}^{-\beta}$ (in units of $\mathrm{erg}\ \mathrm{cm^{-2}}\ \mathrm{s^{-1}}$), where
$\epsilon_{\gamma}$ is the seed photon energy, $E_{\gamma}$ is the energy of the
observed VHE gamma-ray and $\beta$ is the seed photon spectral index. The value of $\beta$ lies in the range $0 < \beta \le 1$.
However, recently, \cite{Sahu:2020tko,Sahu:2020kce} have shown that this standard
photohadronic scenario is inadequate to explain the EHBL-like behavior of the flaring
events of Mrk 501 and 1ES 1959+650. It was also shown that the seed SSC photon flux in
the kinematically allowed region for the $\Delta$-resonance production from $p\gamma$
interaction is no more a single power-law due to the EHBL-like behavior of the flaring
events. So, the standard photohadronic  model was extended to include two different zones for the SSC flux. The zone-1 (low energy region) has the conventional behavior having $0 < \beta \le 1$. But in zone-2 (high energy region) the value of $\beta$ is different and in the range
$ 1 < \beta  \le 1.5$. However, the spectral index $\alpha$ of Fermi accelerated high energy protons in both the zones remain the same.

During the EHBL-like flaring events of Mrk 501 and 1ES 1959+650, the positions of the synchrotron peak and the SSC peak were shifted towards the higher energies. So, \cite{Sahu:2020tko,Sahu:2020kce} assumed that the transition from HBL to EHBL-like behavior is due to the shift of the synchrotron peak and the SSC peak in the SED towards higher energies.
As a result of this shift, jet parameters such as bulk Lorentz factor, blob size, and
magnetic field may be different from HBL, but the particle acceleration and emission
mechanisms remain the same. Therefore, the photohadronic process should still be the
dominant process in this changed environment. It is also argued that these EHBL-like
events are transient in nature and may repeat in the future.

During a MW campaign in 2010, flaring in GeV-TeV energy was observed in Mrk 421. The VHE activity started decreasing from a high flux
to a normal flux during the period from March 10 to March 22, i.e. for 13 consecutive days \citep{Aleksi__2015}. It was
observed that the positions of both the synchrotron peak and the SSC peak of Mrk 421 were shifted towards higher energy values consistent with the EHBL classification of the object, were similar to the
flaring of Mrk 501 and 1ES1959+650 discussed above. This shift in the peak positions towards higher energies is explained using electrostatic acceleration mechanism \citep{Zheng:2020nsk}. Also
recently it has been shown that these VHE spectra of Mrk 421 can be explained better in the
photohadronic model than the one and the two-zone leptonic models \citep{deLeon:2021lof}.  As
the two-zone photohadronic model has successfully explained the EHBL-like behavior of
Mrk 501 and 1ES1959+650, we use this model once again to explain the
EHBL-like behavior of Mrk 421 observed during March 2010.

\section{Flaring of Mrk 421}

Mrk 421 is the nearest ($z=0.031$) well known HBL with a central black hole mass
$M_H\approx (2-9)\times 10^8 M_{\odot}$ and is the first extragalactic source to be
detected in multi-TeV $\gamma$-rays \citep{Punch:1992xw}. It is also one of the fastest
varying $\gamma$-ray sources. Since its discovery, the object has been studied
intensively through dedicated MW observations and several major flares have been
observed \citep{Amenomori:2003gj,Fossati:2007sj,Cui:2004wi,Abeysekara:2016qwu,Mastichiadis:2013aga,Beck:2020gdi}. Its average MW SED has been modelled with both leptonic \citep{Abdo:2011zz,Banerjee:2019oyp} and hadronic \citep{Cerruti:2014iwa,Zech:2017lma} models. During April-May 1994, a TeV/X-ray flare was reported and on May 14/15 a correlation between the TeV emission and X-ray flare was observed. Also, during a
multiwavelength campaign in 2004, large flares in X-ray by Rossi X-ray Timing Explorer (RXTE) and TeV flare by  Whipple 10 m telescope were observed \citep{Blazejowski:2005ih}. It is important to note
that during this period, the TeV flare had no coincident low energy counterparts, and
most significantly,
the X-ray flux reached its peak 1.5 days before the TeV flux during this outburst.
A remarkable similarity between this VHE flare and the orphan TeV flare in 1ES 1959+650
of 2002 \citep{Krawczynski:2003fq} and a similar variation pattern in their X-ray
emission was established. Exceptionally long and dense monitoring in MW of Mrk 421 has
been undertaken since 2009 to understand the temporal evolution of the SED. Mrk 421 was
in active state during a MW campaign in 2010 and an extraordinary flare $\sim 27$ Crab
Unit (c.u.) above 1 TeV was first detected by VERITAS telescopes on February 16 and a
follow-up observation was undertaken by the HESS telescopes \citep{Tluczykont:2011gs}
for four subsequent nights. So far, this is the  highest flare ever observed from Mrk
421. Also, flaring in VHE gamma-ray was observed for 13 consecutive days from March 10
(MJD 55265) to March 22 (MJD55277). Again, during a 6-month long multi-instrument
campaign \citep{Banerjee:2019oyp},
a large VHE flare about 16 times larger than the usual one, was observed by MAGIC
telescopes on 25th of April 2014 and soon it was followed-up by XMM-Newton and VERITAS.
Several of the above multi-TeV flaring events from Mrk 421 have been explained extremely well by the standard photohadronic scenario \citep{Sahu:2015tua,Sahu:2018gik,Sahu:2019kfd}.

\section{Two-zone photohadronic model}
The photohadronic model is very successful in explaining the multi-TeV flaring events from many HBLs
\citep{Sahu:2019lwj,Sahu:2019kfd}. The epoch of VHE flaring is explained by assuming
the formation of a double jet structure along the same axis: an inner jet of size $R'_f$
and photon density $n'_{\gamma,f}$
buried under an outer jet of size $R'_b$ and
photon density $n'_{\gamma}$ where $R'_f<R'_b$ and $n'_{\gamma,f}>n'_{\gamma}$. Here
the notation $'$ implies comoving frame. The internal and the external jets are moving
with almost the same bulk Lorentz factor $\Gamma_{in}\simeq \Gamma_{ext}\simeq \Gamma$.
Also, their common Doppler factor is $\mathcal{D}$ (for blazars $\Gamma \simeq
\mathcal{D}$). The geometrical structure of the jet and the detailed description of the photohadronic model is given in Fig. 1 of \cite{Sahu:2019lwj}.

In the inner jet region the Fermi accelerated protons with a differential power-law spectrum
$dN/dE_p\propto E^{-\alpha}_p$ ($\alpha \ge 2$)
interact with the background SSC seed
photons to produce the $\Delta$-resonance and its subsequent decay gives $\gamma$-rays and neutrinos via the decay of neutral pion and charged pion respectively. As the inner region is hidden, there is no way to directly estimate the photon
density there. Thus we assume a simple scaling behavior of the photon densities in the inner and the outer jet regions satisfying the criteria
\begin{equation}
\frac{n'_{\gamma,f}(\epsilon_{\gamma,1})}{n'_{\gamma,f}(\epsilon_{\gamma,2})} \simeq\frac{n'_{\gamma}(\epsilon_{\gamma,1})}{n'_{\gamma}(\epsilon_{\gamma,2})}.
\label{eq:scal}
\end{equation}
The above equation implies that the ratio of photon densities at two different background energies $\epsilon_{\gamma,1}$ and $\epsilon_{\gamma,2}$ in the inner and the outer jet regions are almost the same. By using the relation in
Eq. (\ref{eq:scal}), we can express the inner photon density in terms of the observed flux. The kinematical condition \citep{Sahu:2019lwj} to produce the $\Delta$-resonance is,
\begin{equation}
E_p \epsilon_\gamma=\frac{0.32\ \Gamma{\mathcal D}}{(1+z)^{2}}\ \mathrm{GeV^2},
\label{eq:kinproton}
\end{equation}
where $z$ is the redshift of the object. The observed VHE $\gamma$-ray energy
$E_{\gamma}$ and the proton energy $E_p$ are related through $E_{\gamma}\simeq
0.1\,E_p$. In previous works by \cite{Sahu:2019lwj} it is explicitly shown that for
the production of $\Delta$-resonance, the value of $\epsilon_{\gamma}$ always lies in
the lower tail region of the SSC band and for HBLs the flux in this region is a perfect power-law.

The VHE $\gamma$-rays en route to Earth undergo
energy dependent attenuation  by the extragalactic background light (EBL)
through electron-positron pair production and this pair-production process
not only attenuates the absolute flux but also significantly changes the
shape of the VHE spectra. The measurement of the EBL is very
difficult due to the uncertainties in the contribution of
zodiacal light \citep{Hauser:2001xs,Chary:2010dc} and consequently
galaxy counts result in a lower limit since the number of unresolved
sources are unknown \citep{Pozzetti:2000ew}. Several approaches with different
degrees of complexity have been developed to
calculate the EBL density as a function of energy for different
redshifts \citep{Dominguez:2010bv,Inoue:2012bk,Franceschini:2008tp,Stecker:2005qs,Stecker:2016fsg}. For the present analysis we use the EBL model of  \cite{Franceschini:2008tp}, which is consistent with a EBL peak flux of $\sim 15\ \mathrm{nW}\mathrm{n^{-2}}\ \mathrm{sr^{-1}}$ at $1.4\, \mathrm{\mu m}$ and is also widely used in other leptonic and hybrid models \citep{Aleksic:2014usa,Ahnen:2016hsf}. Taking into account the EBL contribution the observed
VHE flux from the sources can be given as \citep{Sahu:2019lwj}
\begin{equation}
F_{\gamma,obs}(E_{\gamma})=A_{\gamma} \Phi_{SSC}(\epsilon_{\gamma}) \left ( \frac{E_\gamma}{TeV} \right )^{-\alpha+3}\,e^{-\tau_{\gamma\gamma}(E_\gamma,z)},
\label{eq:fluxphi}
\end{equation}
where $A_{\gamma}$ is a dimensionless constant and $\tau_{\gamma\gamma}$ is the optical depth for the $\gamma\gamma\rightarrow e^+e^-$
process which depends on $E_{\gamma}$ and the redshift.
In the standard photohadronic model, $\Phi_{SSC}\propto E^{-\beta}_{\gamma}$ and putting this in
Eq. (\ref{eq:fluxphi}) the flux can be written as

\begin{equation}
F_{\gamma,obs}(E_{\gamma})=F_0 \left ( \frac{E_\gamma}{TeV} \right )^{-\delta+3}\,e^{-\tau_{\gamma\gamma}(E_\gamma,z)}=F_{\gamma, in}(E_{\gamma})\, e^{-\tau_{\gamma\gamma}(E_\gamma,z)}.
\label{eq:fluxgeneral}
\end{equation}
The normalization constant $F_0$ is deduced from the observed VHE spectrum.
Thus the  spectral index $\delta=\alpha+\beta$ is the only free parameter in this model and
$F_{\gamma,in}$ is the intrinsic VHE $\gamma$-ray flux.
In previous studies \citep{Sahu:2020tko,Sahu:2020kce} it was
shown that the EHBL-like behavior of the flaring epoch can not be explained by
the canonical photohadronic model. In order to remedy this defect the tail region of the SSC band which is
responsible for the production of the $\Delta$-resonance is assumed to have two different
power-laws with spectral indices $\beta_1$ (lower part of the tail region) and $\beta_2$ (upper part of the tail region).
Thus the photon flux is expressed as
\begin{equation}
\Phi_{SSC}\propto
 \left\{
\begin{array}{cr}
E^{-\beta_1}_{\gamma}
, & \quad
 \mathrm{100\, GeV\, \lesssim E_{\gamma} \lesssim E^{intd}_{\gamma}}
\\ E^{-\beta_2}_{\gamma} ,
& \quad   \mathrm{E_{\gamma}\gtrsim E^{intd}_{\gamma}}
\\
\end{array} \right. .
\label{eq:sscflux}
\end{equation}
Here $\beta_1\neq \beta_2$ and $E^{intd}_{\gamma}$ is an energy scale around which the transition takes place and its value can be fixed from the individual flaring spectrum.
We suppose that the change in the value of the spectral index from $\beta_1$ to  $\beta_2$ around $E^{intd}_{\gamma}$ could be due to the effect of magnetic field on the energy loss process of the high energy electrons to the SSC photons or it could be some new physics. However, we do not dwell on the transition between these two zones as the exact physical mechanism behind it is unclear at this stage. The value of $E^{intd}_{\gamma}$ corresponds to a value of $\epsilon_{\gamma}$ in the SSC band. This separates $\Phi_{SSC}$ into two distinct regions. By inserting Eq. (\ref{eq:sscflux}) into Eq. (\ref{eq:fluxphi}), the observed VHE spectrum can be expressed as
\begin{equation}
F_{\gamma, obs}=
e^{-\tau_{\gamma\gamma}}\times
\begin{cases}
 F_1 \, \left ( \frac{E_{\gamma}}{TeV} \right )^{-\delta_1+3}
, & \quad
\mathrm{100\, GeV\, \lesssim E_{\gamma} \lesssim E^{intd}_{\gamma}}\,\,\,\,\, (\text{zone-1}) \\
F_2 \, \left ( \frac{E_{\gamma}}{TeV} \right )^{-\delta_2+3},
& \quad \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \mathrm{E_{\gamma}\gtrsim E^{intd}_{\gamma}}\,\,\,\,\, (\text{zone-2})
 \end{cases}.
\label{eq:flux}
\end{equation}
In this two-zone photohadronic model $F_1$ and $F_2$ are the normalization constants and the spectral indices $\delta_i=\alpha+\beta_i$ ($i=1,2$) are the free parameters to be adjusted by fitting the model to the observed VHE spectrum of the EHBL. The proton spectral index is $\alpha\ge 2.0$, and we use the generally accepted value of $\alpha=2.0$ \citep{Dermer:1993cz}. This automatically constrains the value of $\beta_i$ for a given $\delta_i$. However, it is important to mention that the value of $\delta_1$ (zone-1) allows one to characterize the flaring state for a given observation without depending explicitly on the modeling of the SSC region. In other words, in the photohadronic scenario, the low energy region of the VHE spectrum decides the flaring state of the outburst as the value of $\delta_1$ always lies in the range $2.5\le \ \delta \le 3.0$ \citep{Sahu:2019kfd}.

To fit the observed data for each VHE spectrum, we use a standard optimization of the parameter space through the mean squared error (MSE) with respect to the central values of the data points, and choose the best pair of ($F_{0}, \delta$) by maximizing the Pearson chi-squared ($\chi^2$) statistic and calculating the corresponding statistical significance.

The observed VHE spectrum of each day is fitted in the following manner. For the low energy regime of the spectrum we take $2.5 \le \delta_1 \le 3.0$ and fit it with different values of $F_1$. Similarly in the high energy regime  we take $3.1 \le \delta_2 \le 3.5$ and fit the spectrum for different values of $F_2$. From both the zones we select the best fits by calculating their statistical significance. The crossing point of these two curves gives the value of  $E^{intd}_{\gamma}$. Also, we have to make sure that these two curves intersect only once.

\begin{figure}
\begin{subfigure
\centering
\includegraphics[width=.5\linewidth]{Figure1_MJD55265.pdf}
\end{subfigure}
\begin{subfigure
  \centering
  \includegraphics[width=.5\linewidth]{Figure2_MJD55266.pdf}
\end{subfigure}%
\caption{The VHE flaring of MJD 55265 observed by MAGIC and VERITAS telescopes are shown in (a). Also the VHE event of MJD 55266 observed by MAGIC is shown in (b). Along with the two-zone photohadronic fit (blue for zone-1 and red for zone-2) the one-zone SSC fit for both days are shown for comparison.
In all the figures the normalization constants $F_i$ ($i=1,2$) are defined in units of $\mathrm{erg\, cm^{-2} \, s^{-1}}$. The fitted photohadronic model curves in zone-1 and zone-2 are extended by dashed line of their respective color beyond the transition energy just to show their likely behavior.}
\label{fig:figure1}
\end{figure}
\begin{figure}
\begin{subfigure
\centering
  \includegraphics[width=.5\linewidth]{Figure3_MJD55267.pdf}
\end{subfigure}
\begin{subfigure
\centering
  \includegraphics[width=.5\linewidth]{Figure4_MJD55268.pdf}
\end{subfigure}%
\caption{The VHE flaring data of MJD 55267 observed by VERITAS (a) and of MJD 55268 by both MAGIC and VERITAS (b) are fitted with two-zone photohadronc mdoel which are shown along with the one-zone SSC fit for comparison.}
\label{fig:figure2}
\end{figure}
\begin{figure}
\begin{subfigure
\centering
  \includegraphics[width=.5\linewidth]{Figure5_MJD55269.pdf}
\end{subfigure}
\begin{subfigure
\centering
  \includegraphics[width=.5\linewidth]{Figure6_MJD55270.pdf}

\end{subfigure}%
\caption{The VHE spectra of MJD 55269 (a) and of MJD 55270 (b) observed by both MAGIC and VERITAS telescopes are fitted with two-zone photohadronic model. We have also shown the one-zone SSC fit for comparison.}
\label{fig:figure3}
\end{figure}
\begin{figure}
\begin{subfigure
\centering
  \includegraphics[width=.5\linewidth]{Figure7_MJD55271.pdf}
\end{subfigure}
\begin{subfigure
\centering
  \includegraphics[width=.5\linewidth]{Figure8_MJD55272.pdf}
\end{subfigure}%
\caption{The VHE spectra of MJD 55271 (a) and of MJD 55272 (b) observed by VERITAS telescopes are fitted with two-zone photohadronic model and shown along with their respective one-zone SSC fit for comparison.}
\label{fig:figure4}
\end{figure}
\begin{figure}
\begin{subfigure
\centering
  \includegraphics[width=.5\linewidth]{Figure9_MJD55273.pdf}
\end{subfigure}
\begin{subfigure
\centering
  \includegraphics[width=.5\linewidth]{Figure10_MJD55274.pdf}
\end{subfigure}%
\begin{subfigure
\centering
  \includegraphics[width=.5\linewidth]{Figure11_MJD55276.pdf}
\end{subfigure}%
\begin{subfigure
\centering
  \includegraphics[width=.5\linewidth]{Figure12_MJD55277.pdf}
\end{subfigure}
\caption{The VHE flaring data of MJD 55273 (a), MJD 55274 (b), MJD 55276 (c), and MJD 55277 (d) observed by the Cherenkov telescopes MAGIC and VERITAS are fitted with two-zone photohadronic model. For comparison, we have also plotted the one-zone SSC fits.
}
\label{fig:figure5}
\end{figure}

\section{Results}

During the MW campaign period of 2010, flaring in VHE gamma-rays was observed for 13 consecutive days from March 10 (MJD 55265) to March 22 (MJD55277). Initially the VHE flux was high ($\sim 2$ C.U.) and slowly the flux decreased to a typical value of $\sim 0.5$ C.U. towards the end of the flaring epoch.
In the VHE energy band, MAGIC, VERITAS and Whipple telescope systems observed the flaring events. The MAGIC telescopes made 11 observations and the exposure time for each period was about 10 to 80 min, amounting to a total of 4.7 h of good quality data. However, due to bad weather, data recorded on MJD 55272 and MJD 55275 were removed from the MW observations.

The VERITAS telescopes observed the flaring of Mrk 421 on MJD 55260, 55265, 55267-55274 with a 10 min run time each day. The observations were
performed at zenith angles $18^{\circ}-23^{\circ}$ to benefit from the lowest possible energy threshold of the events. Also, the Whipple 10 m telescope performed ten observations in ON/OFF and TRK (tracking) modes \citep{Pichel:2009gn}, lasting from one to six hours each on MJD 55267-55271 and MJD 55273-55277.

The most important aspect of these 13 days observation is that the synchrotron peak and the SSC peaks were displaced towards higher energies. Particularly, the synchrotron peak frequency was above $10^{17}$Hz and according to the classification scheme of the BL Lac objects, it belongs to EHBL category even though Mrk 421 is well known as an HBL. As discussed in the introduction, similar behavior has also been observed from Mrk 501 and 1ES 1959+650 which are well known HBLs, and there spectra are explained very well with the two-zone photohadronic scenario \citep{Sahu:2020tko,Sahu:2020kce}. Thus using the two-zone photohadronic model, the 13 days VHE flaring events observed by MAGIC and VERITAS are analysed here. For comparison, along with our results the one-zone SSC fits are also shown \citep{Aleksi__2015}. In this analysis, we ignore the statistical significance of a fit when only three or less observed data points are available.

The VHE spectrum of MJD 55265 (10 March 2010) observed by both MAGIC and VERITAS is in the energy range $0.1\, TeV\le E_{\gamma}\le 4.3 \, TeV$ and the synchrotron peak is at $\epsilon_{\gamma}\sim 4.6\times 10^{17}$ Hz which clearly shows its EHBL behavior. Using the two-zone photohadronic model, we fit the low-energy part of the spectrum (zone-1) (for $E^{intd}_{\gamma} \lesssim 1$ TeV) very well with $\delta_1=2.7$ and its corresponding flux normalization is $F_1=2.1\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$, with a statistical significance of 98\% ($\chi^2=0.0008$).
Also the best fit to the high energy part of the VHE spectrum (zone-2) is obtained for $F_2=2.0\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$ and $\delta_2=3.2$, with a statistical significance of 96\% ($\chi^2=0.0366$). These are shown in Fig. \ref{fig:figure1}(a) along with the one-zone SSC model of \citep{Aleksi__2015} for comparison.
It can be seen that our model fits better the data than the one-zone SSC model. The zone-1 and zone-2 have their respective seed SSC flux spectral indices $\beta_1=0.7$ and $\beta_2=1.2$ respectively. The VERITAS data behave differently in its highest energy when the flux has suddenly gone up and is non compatible with the two-zone model. However, previously it was shown that this can be fitted well with conventional photohadronic model with $F_0=2.1\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$ and $\delta=2.6$ corresponding to a very high emission state \citep{Sahu:2019kfd}. As the MAGIC and VERITAS observation windows are seven hours apart, it is possible that intra-night variability might be responsible for this difference.

On MJD 55266 the flaring event was observed by MAGIC telescopes only and the spectrum is in the energy range $0.1\, TeV\le E_{\gamma}\le 4.4$ TeV. The bow shape spectrum on this day (MJD 55266) can not be fitted by a single power-law. The zone-1 ($E^{intd}_{\gamma} \lesssim 1$ TeV) is fitted with $\delta_1=2.7$ and $F_1=2.5\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$ having a statistical significance of 99\% ($\chi^2=0.0189$). Similarly, the best fit to zone-2 is obtained for $F_2=2.45\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$ and $\delta_2=3.5$, with a statistical significance of 82\% ($\chi^2=0.0008$) and on this day the synchrotron peak is at $\epsilon_{\gamma}\sim 5.0\times 10^{17}$ Hz. Our fit and the one-zone SSC results are compared in Fig. \ref{fig:figure1}(b). It can be seen that the leptonic model is unable to fit the data in the low energy regime.

The only data available on MJD 55267 are from VERITAS and the spectrum is very well fitted with the two-zone photohadronic scenario with the transition energy $E^{intd}_{\gamma} \simeq 1$ TeV, which is shown in  Fig. \ref{fig:figure2}(a). The best fitted values for zone-1 are $F_1=2.5\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$ and $\delta_1=2.7$ (statistical significance of 83\%, $\chi^2=0.0538$) while for zone-2 we have $F_2=2.5\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$, and $\delta_2=3.3$.
The SSC model of \cite{Aleksi__2015} predicts a VHE spectrum which falls significantly faster than the observed spectrum.

Both MAGIC and VERITAS observed the flaring event of MJD 55268 and the observed spectra are in the energy range $0.18\, TeV\le E_{\gamma}\le 2.1$ TeV. The MAGIC spectrum is well fitted in zone-1 with the spectral index $\delta_1=2.7$ and the normalization factor $F_1=1.3\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$
and here $E^{intd}_{\gamma} \lesssim 0.8$ TeV .
There are not many data points in zone-2, so we fitted this region with $\delta_2=3.4, 3.5$
and both values fit very well to the observed data. However, above 2 TeV there is a minor difference between these fits. Also our fit to the spectrum differs substantially from the one-zone SSC fit as shown in Fig. \ref{fig:figure2}(b).

The flaring event of MJD 55269 is very similar to the one observed on MJD 55268 and a good fit to the spectrum is achieved with $F_1=2.0\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$  and $\delta_1=2.7$ for zone-1 (statistical significance of 96\%, $\chi^2=0.0790$), while $F_2=1.8\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$, $\delta_2=3.2$ are for zone-2
as shown in Fig. \ref{fig:figure3}(a).

On MJD 55270 the spectrum in the low energy regime i.e. zone-1 is flatter than previous days which implies that the flaring event is moving towards low emission state.
However, because of less number of data points we fitted this region with $\delta_1=2.7, 2.8$
and both fit very well with very similar normalization constants ($F_1=1.32\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$  and $1.16\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$).
On the other hand, zone-2 is fitted well by $F_2=8.7\times 10^{-11}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$, $\delta_2=3.3$ (statistical significance of 87\%, $\chi^2=0.0339$) as shown in Fig. \ref{fig:figure3}(b). One-zone SSC model also fits well to the spectrum. In the high energy regime, however, it falls faster than our fit.

The flaring event of MJD 55271 was only observed by VERITAS in the energy range $0.25\, TeV\le E_{\gamma}\le 2.6$ TeV. Below the transition energy $E^{intd}_{\gamma} \lesssim 0.9$ TeV, the zone-1, we can fit well the spectrum with $F_1=(1.7-1.62)\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$  and $\delta_1=(2.8-2.9)$.
The zone-2 is fitted with $F_2\simeq 1.55\times 10^{-10}\ \mathrm{erg\,cm^{-2}\,s^{-1}}$ and  $\delta_2=3.3$.
The one-zone SSC  model fit is also shown in Fig. \ref{fig:figure4}(a) to compare with our two-zone photohadronic model fit.

For the observations from MJD 55272 to MJD 55277, with the exception MJD 55275 when there was no observation, the value of  $E^{intd}_{\gamma}$ is very low and in a narrow range of  $0.25\, TeV \lesssim E^{intd}_{\gamma}\lesssim 0.4$ TeV. Also the spectrum for each individual day is flat which is best fitted with either $\delta_1=2.9$ or $3.0$ and the statistical significance for these fits ranges from 85\% to 95\%. This value of $\delta_1$ signifies that the flaring of Mrk 421 has reached the low emission state and at the same time the seed photon flux in the SSC band has reached the maximum value of $\beta_1\simeq 1.0$. Also the maximum photon flux in zone-1 has substantially decreased compared to the first day of the observation. The zone-2 for these five days are also fitted very well with $3.2\lesssim \delta_2 \lesssim 3.5$ which correspond to the seed photons in the low energy tail region of the SSC band and their fluxes having spectral indices in the range $1.2\lesssim \beta_2\lesssim 1.5$. The flaring events of all these days along with the one-zone SSC fits are shown in Fig. \ref{fig:figure4}(b) and Figs. \ref{fig:figure5}(a), (b), (c) and (d) respectively.

Initially the VHE flaring was in high emission state and slowly it went down to low emission state by the end of the flaring period. This transition was accompanied by a gradual decrease in the transition energy $E^{intd}_{\gamma}$ from $\sim 1$ TeV to $\sim 0.25$ TeV and consequently zone-2 became wider and spread into zone-1 making the latter region narrower. We recall that the transition in the VHE spectrum took place due to the change in the seed SSC photon flux. So,
physically what it means is that initially the photon flux in the tail region of the SSC band which is responsible for the $\Delta$-resonance production had two distinct energy dependencies as shown in Eq. (\ref{eq:sscflux}) (lower part with $\epsilon^{\beta_1}_{\gamma}$ and upper part with $\epsilon^{\beta_2}_{\gamma}$). Slowly the lower part drifted into the upper part by squeezing the former region into a narrow strip.

On MJD 55266, the maximum observed $\gamma$-ray energy was $E_{\gamma}\simeq 4.4$ TeV which is produced from the interaction of Fermi accelerated proton of energy $E_p\simeq 44$ TeV in the jet with the photons in the lower part of the SSC spectrum. In the two-zone SSC model of \cite{Aleksi__2015} for this day the SSC spectrum starts around $\epsilon_{\gamma}\sim 7\times 10^{21}$ Hz (Figure 8(b) of \cite{Aleksi__2015}). For this value of $\epsilon_{\gamma}$, the minimum value of the bulk Lorentz factor estimated in our model is $\Gamma\simeq 21$. Similar or smaller values for the bulk Lorentz factor are obtained for other days. On the other hand, in the one-zone SSC model (Figure 8(a) of \cite{Aleksi__2015}), the SSC spectrum starts above $2\times 10^{21}$ Hz and for this value of $\epsilon_{\gamma}$ we get minimum
$\Gamma\simeq 35$ which is too high.

\section{Discussion}

Mrk 421 is the nearest HBL and is the first extragalactic source to be observed in VHE. Since then it has undergone several episodes of multi-TeV flaring. Mrk 421 is also the first source observed in multiwavelength for exceptionally long period of time with dense monitoring. The VHE flaring epochs of this source during all its previous observations were consistently found to be in the HBL category which corresponds to the synchrotron peak in the frequency range $10^{15}-10^{16}$ Hz. However, during the 13 consecutive days observations undertaken during March 2010 the source was in the state of high activity with a changed spectral behavior and the synchrotron as well as the second peak shifted towards higher energies. The synchrotron peak shifted above $10^{17}$ Hz, implying that the blazar has undergone a shift from HBL to EHBL. Similar situations are observed in HBLs Mrk 501 and 1ES 1959+650, which were explained very well by the two-zone photohadronic scenario. For Mrk 421, we have used the same two-zone photohadronic model to interpret the 13 consecutive days of flaring in VHE and we find that the individual flaring events are explained very well by our model. This shows once more the success of two-zone photohadronic model in explaining the
extreme HBL nature of the VHE flaring events. We note that when the flaring in VHE was first observed on March 10, 2010 (MJD56265), the source was in high emission state and it continued to be in high emission state for about a week. Then it started degrading to low emission state when the flux was low and also the spectrum in zone-1 was almost flat. During the first week of the observation, the transition energy $E^{intd}_{\gamma}\simeq 1$ TeV when the flaring in the zone-1 was in high emission state. However, from MJD55272, the transition energy $E^{intd}_{\gamma}$ started decreasing and reached a minimum value of $0.25$ TeV on  MJD55274 and for rest of the days the zone-1 was in low emission state. The shift in the $E^{intd}_{\gamma}$ from $\sim 1$ TeV to $\lesssim 0.25$ TeV shows that towards the end of the VHE flaring, the zone-2 had dominated almost the whole spectrum, implying that the seed photon flux in SSC region has changed from a flatter behavior ($\beta_1\simeq 0.7$) to a steeper ($\beta_2\simeq 1.5$) one. This change in spectral behavior, in principle, can be detected by simultaneous observation of the SED in the low energy tail region of the SSC band and in the VHE band during the EHBL-like outburst of the source. We have also estimated the minimum bulk Lorentz factor, obtaining a reasonable value of $\Gamma\approx 21$. Further analysis of EHBL and EHBL-like transient behaviors will help solidify the evidence for the two-zone photohadronic origin of these phenomena.

We are thankful to Luis H. Castañeda Hernandez for many useful
discussions. The work of S.S. is
partially supported by DGAPA-UNAM (Mexico) Project No. IN103019. Partial support from CSU-Long Beach is gratefully acknowledged.

\bibliographystyle{aasjournal}


\section{INTRODUCTION}

\vspace{1mm}

Fuel efficiency is one of the most significant factors for vehicles.
In fact, according to the U.S. Energy Information Administration (EIA) \cite{US_EIA}, nearly 143 billion gallons of motor gasoline were consumed in 2018 in the United States.
Also, $CO_2$ emissions due to motor gasoline and diesel fuel consumption comprise over $30 \%$ of total U.S. energy-related $CO_2$ emissions.

\vspace{1mm}

Connected and Automated Vehicles (CAVs) have immense potential to improve fuel economy.
Even though vehicle automation alone is unlikely to have significant impacts on fuel consumption \cite{wadud2016help}, CAVs can effectively improve fuel efficiency by using intersection connectivity \cite{aoki2017configurable, aoki2019DSIP}, platooning \cite{tsugawa2016review}, eco-driving strategies \cite{peter2020ecodriving} and other novel applications \cite{brown2014analysis}.
With using the peripheral technologies, such as vehicular communications and/or edge computing, CAVs are capable of improving energy efficiency by incorporating a variety of information.

\vspace{1mm}

In particular, Eco-Lane Selection and Eco-Cruise on highways/motorways are promising eco-driving strategies, because there are generally fewer traffic controllers and the vehicles keep moving in general.
In fact, human drivers change the lane when there is a slower lead vehicle and when the neighboring lane has sufficient space, in order to shorten their travel times and save fuel.
The Eco-Cruise maneuver adjusts vehicle speed and acceleration based on traffic, speed limit changes, road features and navigational maneuvers \cite{peter2020ecodriving}.
Eco-Lane Selection determines the lane to drive in and the timing to change the lane, accounting for road structures and traffic.
In these eco-autonomous driving strategies, CAVs incorporate static/dynamic data sources, including the on-board sensors, V2X (Vehicle-to-Everything) communications, pre-loaded vehicle-related data, and map database.
Although V2V (Vehicle-to-Vehicle) communications are not necessary to save fuel, cooperation by V2V communications is at least beneficial for road safety \cite{lehmann2018generic, hess2019negotiation}.

\begin{figure}[!b]
\centering
\includegraphics[width=6.75cm]{pic/MECCconcept.png}
\caption{\textit{MultiCruise} for Energy-Efficient Lane Selection Strategy.}
\label{fig:EcoLaneselection}
\end{figure}

\vspace{1mm}

In this paper, we present a cooperative and energy-efficient lane-selection strategy named {\it MultiCruise} in which each CAV selects the optimal lane and optimal target vehicle to follow among multiple candidates by working with the Eco-Cruise component.
In MultiCruise, as shown in Figure \ref{fig:EcoLaneselection}, the Eco-Cruise component \cite{peter2020ecodriving} calculates multiple sets of the driving parameters of the ego vehicle, and then, the Eco-Lane Selection component determines the lane to drive for saving fuel consumption while ensuring road safety.
MultiCruise accounts for multiple data sources, such as the road curvature, grade, and the surrounding vehicles' velocities and accelerations.
MultiCruise calculates two costs to determine the goal behavior: (i) the cost for staying in the current lane and (ii) the cost for changing the lane and driving the target lane.
The cost function includes fuel costs and penalties.
The fuel costs are for saving fuel, while penalties are for enabling cooperative and collaborative driving to improve road safety during lane-changing maneuvers.
By using these penalties, MultiCruise is able to avoid selfish back-to-back lane-changing maneuvers, which may affect the road safety for the surrounding traffic.
Since CAVs are safety-critical and life-critical applications, MultiCruise never allows a lane-changing maneuver when road safety can be compromised.
Hence, MultiCruise improves fuel economy while satisfying safety constraints.

The primary contributions presented in this paper are as follows:

\vspace{1mm}
\begin{itemize}
\item We present an Eco-Lane Selection strategy named \textit{MultiCruise} that selects the optimal lane and optimal target vehicle to follow to save fuel consumption.
\item We model and develop a cooperative mechanism for a safe lane-changing maneuver that avoids selfish back-to-back lane-changing maneuvers.
\item We test and verify MultiCruise on a co-simulation platform \cite{aoki2020inforich} that includes autonomous driving software and realistic road networks.
\end{itemize}
\vspace{1mm}

The remainder of this paper is organized as follows.
Section I\hspace{-.1em}I reviews previous work related to our research.
Next, Section I\hspace{-.1em}I\hspace{-.1em}I describes the lane-changing maneuvers focused in this paper.
Section I\hspace{-.1em}V presents our eco-lane selection strategies.
Section V gives the implementation and evaluation of our eco-driving strategies.
Finally, Section V\hspace{-.1em}I presents our summary and future work.

\vspace{0.5mm}

\section{RELATED WORK}

The advent of Connected and Automated Vehicles (CAVs) \cite{urmson2008autonomous} is a great opportunity to improve not only safety but also vehicle fuel economy.
In fact, ARPA-E runs NEXTCAR program (Next-Generation Energy Technologies for Connected and Automated On-Road Vehicles) that aims to leverage the CAV technologies to further improve the energy efficiency \cite{aoki2020inforich, olin2019reducing, oncken2020connected}.
Vehicle automation alone is unlikely to have significant impacts on fuel savings \cite{wadud2016help}, but there are multiple promising technologies for fuel savings, such as vehicle platooning, intersection connectivity, traffic management and eco-autonomous driving strategies.

In particular, eco-autonomous driving strategies can effectively improve fuel efficiency and they include four eco-driving applications: Eco-Approach, Eco-Departure \cite{barth2019ead, shao2020eco}, Eco-Cruise \cite{peter2020ecodriving, jing2016design, firoozi2019safe} and Eco-Lane Selection \cite{aoki2020towards}.
Eco-Approach is defined as the fuel-efficient vehicle/powertrain operation that applies smooth deceleration to bring a vehicle to stop in fuel-friendly fashion.
The Eco-Departure maneuver should conduct a smooth and less-aggressive acceleration profile to reach the target cruise speed. The target speed and the distance to reach the speed typically depend on the traffic preview information that are available from on-board sensors and vehicular communications. Eco-Approach and Eco-Departure are closely related to road intersections and traffic controllers \cite{barth2019ead}, especially signalized intersections.
For example, by using the SPaT (Signal Phase and Timing) information via vehicular communications, each connected vehicle can save fuel.

In addition, Eco-Cruise is the vehicle velocity control strategy to save the fuel consumption.
One of the promising techniques for the cruise control is a Model Predictive Control (MPC), and Jing et al. \cite{jing2016design} used a Finite State Machine (FSM) to narrow down the search space for each MPC step to be used for real-time applications.
Also, Firoozi et al \cite{firoozi2019safe} presented a Eco-Cruise component that incorporates the road grade information and V2V communications.
In \cite{peter2020ecodriving}, Eco-Cruise is designed and developed to determine the driving parameters of the ego vehicle by using the given lane and given target vehicle at the moment. The Eco-Cruise component accounts for road grade, curvature and surrounding traffic, and we use this component for MultiCruise.

Eco-Lane Selection is also a significant application for fuel savings.
FWHA's Next Generation SIMulation (NGSIM) program \cite{NGSIM2006} has defined ``target lane'' concept and also has introduced Freeway Lane Selection (FLS) algorithm, for modeling a lane-changing maneuver on highways and/or motorways.
For connected vehicles, Jin et al \cite{jin2014improving}, Kang et al. \cite{kang2019development} and Tian et al. \cite{tian2018connected} proposed optimal lane-selection algorithms for traffic operations and management on highway segments. Since these works relied on V2I (Vehicle-to-Infrastructure) communications and segment-based operations, the approaches were not scalable. Also, there are multiple on-going studies for lane-changing and merging maneuvers \cite{aoki2018dynamic, aoki2017merging} for autonomous vehicles.
Unlike our work presented in this paper, these works just focused on the throughput and never discussed fuel efficiency.

\section{MANDATORY AND DISCRETIONARY LANE-CHANGING MANEUVERS}

In this section, we present a review of lane-changing maneuvers on highways.
Lane-changing maneuvers have been originally studied for human-driven vehicles \cite{jin2014improving, goswami2007gap, zhao2017analysis}.
These maneuvers can be categorized into two classes \cite{toledo2003modeling}: (i) \textbf{Mandatory lane-changing} and (ii) \textbf{Discretionary lane-changing}.
First, mandatory lane-changing maneuvers are conducted for vehicle navigation.
When a vehicle has to change its lane to drive, for lane reduction or for exiting from highways, the lane-changing maneuver is classified into mandatory one.
For such mandatory lane-changing maneuvers, when the vehicle cannot complete its maneuver until the deadline, it may have to re-calculate its route by using the navigation system.
On the other hand, discretionary lane-changing maneuvers are for secondary purposes, such as energy efficiency, traffic throughput and/or passenger comfort.
A vehicle does not need to complete the discretionary lane-changing maneuver for navigation purposes.
Several researchers have studied the lane-changing maneuvers for autonomous driving \cite{dong2018continuous, krasowski2020safe}, but they do not focus on such features \cite{dong2018continuous, krasowski2020safe}.

\vspace{1mm}

In this paper, we focus on discretionary lane-changing maneuvers, because our primary objective is saving fuel consumption on highways and motorways.
Here, MultiCruise for CAVs is the lane-selection component for saving energy and it is not for vehicle navigation and route planning. Therefore, each CAV only uses MultiCruise to change its lane when it can follow the route provided by its navigation system.

\section{MULTICRUISE: AN ECO-LANE SELECTION STRATEGY}

We now present MultiCruise as an energy-efficient lane-selection strategy.
MultiCruise consists of $2$ steps to determine the lane to drive for saving fuel consumption: (i) \textit{Calculating cost for each potential target} and (ii) \textit{Determining goal behaviors}.
In this section, we first present the cost function, and secondly propose a decision-making policy.
In addition, to show the safety and practicality of MultiCruise, we discuss cooperative and safe mechanisms for MultiCruise.

Each CAV uses MultiCruise by incorporating various data sources, including the on-board sensors, V2X (Vehicle-to-Everything) communications, pre-loaded vehicle-specific lookup tables, and map database.
Since each vehicle determines its behaviors in a distributed manner, our MultiCruise is designed to be used both around human-driven vehicles and around CAVs seamlessly.

\subsection{Cost and Penalty: For Multiple Targets}

Each CAV uses MultiCruise to calculate three Eco-Driving costs to determine whether it should change the driving lane for saving fuel consumption. These three Eco-Driving costs are captured in Eq. (\ref{eq:three_cost}).
Here, we calculate the driving costs on the left neighboring lane $C_{LF}$, current lane $C_{CF}$, and right neighboring lane $C_{RF}$, by using the cruise cost $\Gamma$ introduced and discussed in \cite{peter2020ecodriving} and the penalty $P$ for safety purposes.
The penalty $P$ is presented in Eq. (\ref{eq:penalty_function}).

\begin{eqnarray}
\left\{ \begin{array}{lll}
C_{LF} &= & \Gamma(a_{Ego}, v_{Ego}, v_{LF}, d_{LF}) + P(t_\Delta, v_{Ego})\\
C_{CF} &= & \Gamma(a_{Ego}, v_{Ego}, v_{CF}, d_{CF}) \\
C_{RF} &= & \Gamma(a_{Ego}, v_{Ego}, v_{RF}, d_{RF}) + P(t_\Delta, v_{Ego})\\
\end{array} \right.
\label{eq:three_cost}
\end{eqnarray}

\begin{equation}
P(t_\Delta, v_{Ego}) = s \cdot \frac{v_{Ego}}{t_\Delta}
\label{eq:penalty_function}
\end{equation}

Here, $a_{Ego}$ and $v_{Ego}$ are the ego vehicle's acceleration and speed, respectively.
$d_{LF}$, $d_{CF}$ and $d_{RF}$ represent the inter-vehicle distances to the lead vehicle on the left lane, current lane, and right lane, respectively, as shown in Figure \ref{fig:EcoLaneselection}.
Also, $v_{LF}$, $v_{CF}$ and $v_{RF}$ are the lead vehicles' speeds on the left lane, current lane, and right lane, respectively.
The cruise cost $\Gamma$ \cite{peter2020ecodriving} is calculated in the Eco-Cruise component and it includes (i) fuel cost, (ii) progress cost, and (iii) comfort cost.
The fuel cost is calculated by using the vehicle-specific fuel table, to get the expected fuel consumption from the vehicle velocity and acceleration.
The progress cost is designed to penalize when the ego vehicle is too slow and does not make progress over time.
The comfort cost accounts for vehicle jerk, in order to provide comfortable driving to the human passengers.
A more detailed discussion on jerk is presented in \cite{peter2020ecodriving}.

In addition, the penalty $P$ in Eqs. (\ref{eq:three_cost}) and (\ref{eq:penalty_function}) is designed for avoiding back-to-back and aggressive lane-changing maneuvers.
Here, $t_\Delta$ describes the elapsed time from the last lane-changing maneuver. $s$ is a scaling factor.
By using this mechanism, the penalty becomes very large value right after the lane-changing maneuver is completed, and consequently, MultiCruise does not encourage the lane-changing maneuver at this moment.
The vehicle may be able to reduce fuel consumption by conducting back-to-back lane-changing maneuvers, but this behavior compromises road safety \cite{li2006realistic} and we do not encourage this behavior in MultiCruise.

Each CAV uses these $3$ Eco-Driving costs, $C_{LF}$, $C_{CF}$ and $C_{RF}$, to determine its behavior. We will propose and discuss the decision-making policy next.

\begin{figure}[!t]
\centering
\includegraphics[width=8.55cm]{pic/lanechanging_3.png}
\caption{Space Checker for Lane-Changing Maneuvers.}
\label{fig:Lane_selection_target}
\end{figure}

\subsection{Decision-making Policy}

MultiCruise compares the $3$ Eco-Driving costs and selects one lane to save fuel consumption.
To guarantee road safety under the practical CAV systems, MultiCruise has to satisfy $2$ requirements: (i) MultiCruise does not overturn the decision of the navigation system, and (ii) MultiCruise enables a safe lane-changing maneuver.
For the first requirement, we use a threshold distance $D_{Thr}$ and a distance for the navigation $D_{NAV}$.
$D_{NAV}$ represents the distance to the end of the road segment and/or ramp-off.
When $D_{NAV}$ is smaller than $D_{Thr}$, each vehicle has to consider the mandatory lane-changing maneuver for its navigation purpose.
For the second requirement, after determining the target lane $\theta$, each vehicle has to confirm that the target lane has the sufficient space to complete its lane-changing maneuver by using the sensor-based perception. Since the lane-changing maneuvers may lead vehicle accidents and/or collisions, each vehicle can start the lane-changing maneuver only when safety is guaranteed.

The decision-making algorithm for MultiCruise is captured in Algorithm \ref{algorithm:decition-making}.
Here, we define and use $C_{Thr}$ as a threshold value for the cost, in order to avoid unnecessary lane-changing maneuvers, which may result only in trivial benefits.
Since the lane-changing maneuvers may disturb the entire traffic flow, MultiCruise does not encourage the aggressive and/or unnecessary lane-changing maneuvers.

\begin{algorithm}[b!]
    \If{$D_{NAV} > D_{Thr}$}
    {
      Calculate $C_{LF}$, $C_{CF}$, $C_{RF}$\;
      \If{$C_{CF}>C_{Thr}$}
      {
       target lane $\theta = \mathop{\rm arg~min}\limits_{\alpha \in \{LF, LC, LR\}} C_{\alpha}$\;
      }
    }
\caption{Decision-making algorithm for MultiCruise}
\label{algorithm:decition-making}
\end{algorithm}

\subsection{Safety for Lane-changing Maneuver}

Although the primary objective of MultiCruise is saving fuel consumption, to enhance the road safety for lane-changing maneuvers, MultiCruise includes $2$ mechanisms: (i) Passed time $t_\Delta$ and (ii) Space checker.

First, $t_\Delta$ represents the passed time from the last lane-changing maneuver.
When the value $t_\Delta$ is small, the penalty $P$ becomes a large value and MultiCruise discourages the lane-changing maneuver at this moment.
Therefore, after each vehicle changes its lane to drive, it stays in that lane for a certain time period.
Since back-to-back lane-changing maneuvers disturb the surrounding traffic and may lead to collisions, this cooperative mechanism enhances road safety.

Secondly, when each vehicle uses MultiCruise, it always confirms the space on the target lane, as shown in Figure \ref{fig:Lane_selection_target}, by using its space checker.
In Figure \ref{fig:Lane_selection_target}, the green box represents the area confirmed by the space checker. When there is a moving obstacle within the area, MultiCruise stops the lane-changing maneuver, because it cannot guarantee road safety.
MultiCruise uses the space checker after the target lane $\theta$ is selected by the decision-making algorithm described in Algorithm \ref{algorithm:decition-making}.

By using these mechanisms, MultiCruise improves road safety while enhancing fuel efficiency on highways and/or motorways.

\begin{figure}[!b]
\centering
\includegraphics[width=8.55cm]{pic/Hardwaresetup.png}
\caption{Co-simulation Hardware Setup.}
\label{fig:harware_setup}
\end{figure}

\begin{figure}[!b]
  \begin{center}
    \begin{tabular}{c}

      \begin{minipage}{0.45\hsize}
        \begin{center}
          \includegraphics[width=4.15cm]{pic/vires_moving.png}
	\hspace{7.2cm} (a) Moving Obstacles \\ in Traffic Simulator.

        \label{fig:Synchro_2size}
        \end{center}
      \end{minipage}

      \begin{minipage}{0.55\hsize}
        \begin{center}
          \includegraphics[width=4.15cm]{pic/iread_moving2.png}
	\hspace{7.2cm} (b) Moving Obstacles \\ in iREAD System.

          \label{fig:iread_moving}
        \end{center}
      \end{minipage}\\ \\

    \end{tabular}
    \caption{InfoRich Co-simulation Environment.}
\label{fig:overview_co-simsystems}
  \end{center}
\end{figure}

\begin{figure}[!t]
\centering
\includegraphics[width=7.55cm]{pic/MultiCruise_System2.png}
\caption{MultiCruise Implementation on the Co-simulation Platform \cite{aoki2020inforich}.}
\label{fig:system_implementation}
\end{figure}

\begin{figure}[!t]
  \begin{center}
    \begin{tabular}{c}

      \begin{minipage}{0.5\hsize}
        \begin{center}
          \includegraphics[width=3.9cm]{pic/LongNREL_A_VIRES.png}
          \hspace{6.8cm} (a) \textit{Map A}: Road networks.
          \label{fig:LongNREL_A}
          \end{center}
      \end{minipage}

      \begin{minipage}{0.5\hsize}
        \begin{center}
          \includegraphics[width=3.9cm]{pic/LongNREL_D_VIRES.png}
          \hspace{6.8cm} (b) \textit{Map D}: Road networks.
          \label{fig:LongNREL_B}
          \end{center}
      \end{minipage}
    \end{tabular}
    \caption{Synthetic Drives with a Highway Road Segment.}
    \label{fig:LongTrip_NREL}
  \end{center}
\end{figure}
\section{IMPLEMENTATION AND EVALUATION}

In this section, we present the implementation and evaluation of our MultiCruise in the InfoRich co-simulation platform introduced in \cite{aoki2020inforich} that includes autonomous driving software \cite{bhat2018tools} and a Simulink-based VD\&PT model, as shown in Figure \ref{fig:harware_setup}. In the co-simulation platform, each simulator keeps exchanging the data by using the UDP protocol.
The GUIs for InfoRich co-simulation platform is presented in Figure \ref{fig:overview_co-simsystems}.
The traffic and SPaT information are generated in the traffic simulator as shown in Figure \ref{fig:overview_co-simsystems}-(a), and these information is processed to determine vehicle's behaviors in the autonomous driving software as shown in Figure \ref{fig:overview_co-simsystems}-(b).
Such co-simulation is a promising technique to model complex and dynamic systems in a distributed manner. In the co-simulation system, each subsystem runs without being aware of the entire system.
We evaluate the energy-efficient lane selection strategy in terms of fuel efficiency, and compare against a baseline protocol.
At the same time, we show the feasibility, safety and practicality of the iREAD co-simulation platform \cite{aoki2020inforich} by using it for testing and verifying our eco-driving strategy.

\begin{figure*}[!t]
  \begin{center}
    \begin{tabular}{c}

      \begin{minipage}{0.33\hsize}
        \begin{center}
          \includegraphics[width=5.25cm]{pic/MultiCruise_Distance.jpg}
          \hspace{7.8cm} (a) Travel Distance for Simple Scenario.
          \label{fig:SimpleDistance}
          \end{center}
      \end{minipage}

      \begin{minipage}{0.33\hsize}
        \begin{center}
          \includegraphics[width=4.75cm]{pic/MultiCruise_Yaw.jpg}
          \hspace{6.8cm} (b) Yaw Changes for Simple Scenario
          \label{fig:SimpleYawChanges}
          \end{center}
      \end{minipage}

      \begin{minipage}{0.33\hsize}
        \begin{center}
          \includegraphics[width=4.65cm]{pic/MultiCruise_FuelTotal.jpg}
          \hspace{6.8cm} (c) Consumed Fuel for Simple Scenario
         \label{fig:SimpleFuelTotal}
          \end{center}
      \end{minipage}
    \end{tabular}
    \caption{Performance Evaluations of MultiCruise for a Simple Scenario.}
    \label{fig:SimpleEvaluations}
  \end{center}
\end{figure*}

\begin{figure*}[!t]
\centering
\includegraphics[width=19.0cm]{pic/MultiCruise_DuelRate.jpg}
\caption{Fuel Rates for a Simple Scenario.}
\label{fig:SimpleFuelRate}
\end{figure*}

First, we show the implementation on the eco-autonomous driving system in Figure \ref{fig:system_implementation}.
In the system, the Route Planner component calculates the route to be taken, by using the given sensor information and map information.
Based on the output, Behavioral Planner determines the lane to drive, and this component includes MultiCruise.
In Behavioral Planner, the Eco-Behavior Goals components transmit the information of the left lane, current lane and right lane, and then Eco-Trajectory Generation \cite{peter2020ecodriving} calculates the speed trajectories and cruising costs $\Gamma$ for these three lanes.
Finally, by using these three values, Eco-Behavior Goals determines the lane to drive while accounting for safety requirements.

We first evaluate our eco-driving strategy with a simple scenario in which there is only a straight road segment. In this scenario, we show the details of the vehicle movements, fuel rates, and total fuel consumption through each drive.
In addition, we evaluate MultiCruise with realistic synthetic scenarios discussed in \cite{aoki2020inforich}.
We prepare $6$ different road networks that have a variety of features, such as speed limit, road grade, curvature, and traffic controllers. These synthetic scenarios have a highway road segment, and they have more than $12$ (km) of total distance traveled.

\subsection{Performance Metric}

To evaluate our MultiCruise in terms of fuel efficiency, we use \textit{Consumed Fuel} as a metric.
The Consumed Fuel represents the fuel consumption used in each drive. The primary purpose of MultiCruise is to decrease Consumed Fuel.
The Consumed Fuel is mainly determined by vehicle dynamics and powertrain models, and the co-simulation test platform \cite{aoki2020inforich} enables to evaluate the performance.

\begin{table}[b]
\centering
\caption{Basic Information for Synthetic Cycles.}
  \begin{tabular}{c|c|c|c|c
                          & Overall & Highway       & Speed Limit & Number of\\%        & Ave \#  & Max \# & Min \# & Ave \# & Max \# & Min \#\\
                          & Distance & Distance & on Highway & Intersectons\\ \hline \hlin
    \footnotesize{\textit{Map A}} & 20.31 (km) & 12.71 (km) & 105 (km/h) & 15\\ \hline
    \footnotesize{\textit{Map B}} & 17.33 (km) & 9.86 (km) & 105 (km/h) & 10\\ \hlin
    \footnotesize{\textit{Map C}} & 19.99 (km) & 4.89 (km) & 97 (km/h) & 26\\ \hlin
    \footnotesize{\textit{Map D}} & 15.06 (km) & 10.79 (km) & 105 (km/h) & 8\\ \hlin
    \footnotesize{\textit{Map E}} & 15.59 (km) & 3.72 (km) & 112 (km/h) & 10\\ \hlin
    \footnotesize{\textit{Map F}} & 12.71 (km) & 3.46 (km) & 112 (km/h) & 14\\
  \end{tabular}
\label{table:NREL_highway}
\end{table}

\subsection{Synthetic Scenarios}

To show the feasibility and practicality of MultiCruise, we utilize synthetic drive cycles derived from real-world driving data provided from the National Renewable Energy Laboratory (NREL).
In particular, we select $6$ synthetic drives, \textit{Map A} to \textit{Map F}, that include highways or motorways, because MultiCruise focuses on the discretionary lane-changing maneuvers on highways.
We present $2$ examples of the synthetic cycles, as shown in Figure \ref{fig:LongTrip_NREL}.
Also, in Table \ref{table:NREL_highway}, we present the basic information of the $6$ synthetic cycles, including the cycle distance, the highway distance, the speed limit on the highway segment, and the number of road intersections.
For example, as shown in Table \ref{table:NREL_highway}, \textit{Map A} has more than $20$ (km) for overall distance and includes approximately $12$ (km) along a highway segment, in which there are no traffic controllers and/or intersections.

\subsection{Evaluation with a Simple Scenario}

We first evaluate our MultiCruise with a simple straight-road segment.
The road segment has multiple lanes to drive and has $1,600$ meters.
Also, for this case study, the speed limit is set as $112$ (km/h) that is equal to $70$ (miles/h).
We evaluate the fuel consumption, compared to a baseline \textbf{Lane-following Protocol}.
In the Lane-following Protocol, each vehicle follows the lane to drive for its navigation purpose and does not change its lane for fuel saving and/or traffic throughput.

In addition, for this case study, we have a slower lead vehicle that drives approximately at $70$ (km/h).
In the lane-following protocol, each vehicle simply drives and stays in its current lane even when there is a slower lead vehicle.

We present the traveled distance, yaw angle changes and total fuel consumption in Figures
\ref{fig:SimpleEvaluations}-(a), -(b) and -(c), respectively.
Also, we present the fuel rate in Figure \ref{fig:SimpleFuelRate}.
First, as shown in Figure \ref{fig:SimpleEvaluations}-(a), MultiCruise allows the ego vehicle to reach the destination in a much shorter trip time, because the vehicle can adaptively select the lane and start the lane-changing maneuver when it detects the slower lead vehicle.
On the other hand, the baseline protocol spends a longer trip time because of the slow lead vehicle.
Secondly, we show the yaw changes of the travels with MultiCruise and with the baseline protocol. The vehicle with MultiCruise changes its lane twice to overtake the slower lead vehicle.
On the other hand, the vehicle under the baseline protocol does not change its driving lane, and therefore, the yaw angle is constant.

Overall, MultiCruise saves around $32 \%$ fuel in this case study, as shown in Figure \ref{fig:SimpleEvaluations}-(c). The vehicle with the lane-following protocol has to slow down once it is behind the slower lead vehicle, and consequently, it wastes energy. The vehicle with MultiCruise not only does not need to hit the brakes but also can avoid wasting fuel by dynamically selecting the lane to drive at a higher speed.

\vspace{1mm}
Since this case study is with a simple scenario, we study MultiCruise in more scenarios next.

\begin{figure}[!t]
\centering
\includegraphics[width=8.00cm]{pic/Fuel_LaneChange2.png}
\caption{Performance Evaluation of MultiCruise on the Fuel Consumption for Synthetic Drives.}
\label{fig:Synthetic_Fuel}
\end{figure}

\subsection{Evaluation with Synthetic Scenarios}

We evaluate MultiCruise with $6$ synthetic scenarios, \textit{Map A} to \textit{Map F}.
Each road network contains $1$ highway segment.

We compare MultiCruise to the lane-following protocol under $2$ traffic volumes: (i) Moderate traffic (100 (m/vehicle)) and (ii) Heavy traffic (50 (m/vehicle)). Since MultiCruise focuses on lane-changing maneuvers for energy efficiency, we do not evaluate the performance under the traffic jam and/or congestion.
Under the serious traffic jam, since vehicles move very slowly, the eco-lane selection strategy does not improve fuel efficiency.

As shown in Figure \ref{fig:Synthetic_Fuel}, we present the relative fuel consumption in moderate traffic and in heavy traffic.
We show the relative performance of MultiCruise, compared to the baseline protocol.
As shown in the result, MultiCruise has superior performance and successfully saves Consumed Fuel across all the cases.
In particular, MultiCruise saves fuel consumption by $7.2 \%$ in \textit{Map A} and by $8.5 \%$ in \textit{Map B} with moderate traffic.
In \textit{Map C} and \textit{Map F}, MultiCruise has a small improvement in fuel efficiency, because the highway segment is relatively short and there are few opportunities to change lanes.
One interesting point to note is that MultiCruise saves more energy in moderate traffic than in heavy traffic, other than in \textit{Map D}.
In moderate traffic, there are more spaces and opportunities to conduct lane-changing maneuvers, and MultiCruise contributes to saving fuel consumption.
As discussed and shown in this section, MultiCruise improves fuel efficiency both in a simple scenario and in practical cases.

\section{SUMMARY AND FUTURE WORK}

In this paper, we presented a cooperative and energy-efficient lane-selection strategy named \textit{MultiCruise}, where each Connected and Automated Vehicle (CAV) selects the optimal lane and optimal target vehicle to follow from multiple candidates to save fuel consumption.
We described the Eco-Driving costs to determine whether the vehicle should stay or change the driving lane for fuel standpoint.
Also, we presented safe and cooperative mechanisms for MultiCruise that avoid selfish back-to-back lane-changing maneuvers and avoid vehicle collisions and/or accidents.
MultiCruise incorporates the Eco-Cruise component that adjusts vehicle speed and acceleration for eco-autonomous driving.
We conclude that MultiCruise improves fuel efficiency on highways and motorways by showing both simple case study and realistic simulation with the co-simulation test platform and the autonomous driving software.
In future work, we will study the energy optimization across multiple CAVs, because our work focused on the energy consumption only for the ego vehicle.
In addition, we will develop a Hardware-in-the-Loop (HIL) simulator by integrating with a real vehicle and implement MultiCruise on it.

\vspace{1mm}

\section*{Acknowledgment}
This paper is based upon the work supported by the United States Department of Energy (DOE), ARPA-E NEXTCAR program under award No. DE-AR0000797.
The authors would like to thank Jeff Gonder and Michael O'Keefe at National Renewable Energy Laboratory (NREL) for providing real world driving data and analysis.

\bibliographystyle{ieeetr}


\section{Introduction}

With the booming of deep learning in natural language processing (NLP) field, many pre-trained language models (PLMs) are proposed, such as BERT~\cite{bert}, XLNet~\cite{xlnet}, RoBERTa~\cite{roberta}, ALBERT~\cite{albert}, T5~\cite{t5}, ELECTRA~\cite{electra} and so on, achieving state-of-the-art (SOTA) performance on various tasks. However, these large PLMs are computationally expensive and require a large memory footprint, which makes it difficult to execute them on resource-restricted devices.

To tackle this problem, many works propose task-agnostic BERT distillation~\cite{distilbert, tinybert, minilm} to obtain a general small BERT model that can be fine-tuned directly as the teacher model (e.g., BERT-Base) does. However, the process of task-agnostic BERT distillation is also computationally expensive~\cite{structuredqa}. Because the corpus used in the distillation is large-scale, and each training step is computationally consuming that a forward process of teacher model and a forward-backward process of student model need to be performed.

\begin{figure}[t]
\includegraphics[width=0.48\textwidth]{etd_v8.pdf}
\caption{Overview of our ETD strategy.}
\centering
\label{fig:etd}
\end{figure}

Reusing the parameters of the teacher model to initialize the student model has been proved effective to improve the efficiency of BERT distillation (e.g., BERT-PKD~\cite{bertpkd} and DistilBERT~\cite{distilbert}). However, these works only do the \emph{depth-wise} extraction that simply copies some Transformer layers [22] from the teacher model, which requires the student to keep the same setting of \emph{width dimensions} including hidden/head/FFN dimension, as the teacher. Actually, in many industrial scenarios, we need to build the student models with different widths and depths to meet various latency and memory requirements.

In this paper, we propose \emph{Extract Then Distill} (ETD) as shown in Fig.~\ref{fig:etd}, a flexible and effective method to reuse the teacher's parameters to initialize the student, which firstly allows the student to have a narrower width than the teacher. Specifically, we propose two methods: ETD-Rand and ETD-Impt, which width-wise extract the teacher's parameters randomly and depending on the importance scores respectively. Then we adopt the strategy of uniform layer selection for depth-wise extraction. Finally, we initialize the student with the extracted parameters, then perform the task-agnostic distillation.

Our contributions are two-fold: (1) we propose an effective method ETD, which improves the efficiency of task-agnostic BERT distillation by reusing the teacher parameters to the initialize the student model. The proposed ETD method is flexible and applicable to student models of any size. (2) We demonstrate the effectiveness of ETD on the GLUE benchmark and SQuAD. The extensive experimental results show that ETD can save 70\% of computation cost of the baseline, and when using the same computing resources, ETD outperforms the baselines. ETD is general and can be applied to different existing state-of-the-art distillation methods, such as TinyBERT and MiniLM, to further boost their performance. Moreover, the extraction process of ETD is efficient and brings almost no additional calculations.

\section{Preliminary}
\subsection{Architecture of BERT}
\noindent {\bf Embedding Layer.} Through embedding layer, all tokens in a sentence are mapped to vectors of the \textit{hidden size} $\{\bm{x}_i\}^{|x|}_{i=1}$, where $|x|$ means the number of tokens in this sentence. A special token \textsc{[CLS]} is added at the beginning of the sentence for obtaining the representation of the entire sentence.

\noindent {\bf Stacked Transformer Layers.} Transformer layer dynamically encodes context information into the representation vector of each token through the self-attention mechanism, which tackles the problem of semantic ambiguity to a great extent. We denote the output of embedding layer $[\bm{x}_1;...;\bm{x}_{|x|}]$ as \emph{hidden states} $\bm{H}^0$. The stacked Transformer layers compute the contextual vectors as:
\begin{equation}
    H^l = {\rm Transformer}_l(H^{l-1}), l\in \left[1, L\right]
\end{equation}
where L means the number of Transformer layers.

\subsection{Architecture of Transformer}
Each Transformer layer consists of two sub-modules: the multi-head attention (MHA) and the fully connected feed-forward network (FFN). \emph{Residual connection} and \emph{layer normalization} are employed on top of each sub-module.

\noindent {\bf MHA.} In each layer, Transformer uses multiple self-attention heads to aggregate the output vectors of the previous layer. For the $(l+1)$-th Transformer layer, the output of MHA is computed via:
\begin{equation}
\begin{split}
    &\bm{Q}_i=\bm{H}^l\bm{W}_i^Q,\bm{K}_i=\bm{H}^l\bm{W}_i^K,\bm{V}_i=\bm{H}^l\bm{W}_i^V,\\
    &\bm{A}_i=\frac{\bm{Q}_i\bm{K}_i^T}{\sqrt{d_k}},\\
    &{\rm head}_i = {\rm softmax}(\bm{A}_i)\bm{V}_i,\\
    &{\rm MHA}(\bm{H}^l) = {\rm Concat(head_1, ..., head_{\it a})}\bm{W}^O \\
    &\bm{H}^{\rm MHA} = {\rm LayerNorm}(\bm{H}^l + {\rm MHA}(\bm{H}^l)).
\end{split}
\end{equation}

The previous layer’s output $\bm{H}^l$ is linearly projected to queries, keys and values using parameter matrices $\bm{W}^Q, \bm{W}^K, \bm{W}^V$ respectively. ${\rm head}_i$ indicates the context-aware vector which is obtained by the scaled dot-product of queries and keys in the $i$-th attention head. $a$ represents the number of self-attention heads. $d_k$ is the dimension of each attention head acting as the scaling factor, named as head dimension. $d_k \times a$ is equal to the hidden dimension $d_h$. $\bm{H}^{\rm MHA}$ will be the input of the next sub-module FFN.

\noindent {\bf FFN.} It consists of two linear layers and a GeLU~\cite{gelu} function between them:
\begin{equation}
\begin{split}
    &{\rm FFN}(\bm{H}) = {\rm GeLU}(\bm{H}\bm{W}_1 + \bm{b}_1)\bm{W}_2 + \bm{b}_2,\\
    &\bm{H}^{l+1} = {\rm LayerNorm}(\bm{H}^{\rm MHA} + {\rm FFN}(\bm{H}^{\rm MHA})).
\end{split}
\end{equation}

\subsection{Overview of BERT Distillation}\label{subsec:transdil}

There are two widely-used distillation layer mapping strategies, \emph{uniform} and \emph{last-layer}, adopted by TinyBERT~\cite{tinybert} and MiniLM~\cite{minilm} respectively. We adopt the last-layer strategy and follow other settings of TinyBERT, because we empirically find that the last-layer strategy performs better. The last-layer strategy aims to minimize the mean squared errors between the attention patterns $\bm{A}^S$, output hidden states $\bm{H}^S$ of the student’s last Transformer layer and the counterpart $\bm{A}^T,\bm{H}^T$ of the teacher’s last Transformer layer. The overall objective of the last-layer distillation is formulated as:

\begin{equation}
    \mathcal{L}_{\rm KD} = \frac{1}{a} \sum_{i=1}^{a}{\rm MSE}(\bm{A}_i^S, \bm{A}_i^T) + {\rm MSE}(\bm{H}^S\bm{W}, \bm{H}^T),
\end{equation}
where $\bm{A}_i^S$ and $\bm{A}_i^T$ mean the attention matrix corresponding to the $i$-th head in the last layer of student and teacher respectively. The learnable matrix $\bm{W}$ is introduced to solve the problem of mismatched hidden size between student model and teacher model.
\section{Methodology}
Our proposed ETD strategy is general and can be applied to students of any size. It consists of three steps as shown in Fig.~\ref{fig:etd}: (i) \textbf{Width-wise Extraction} that extracts parameters from the teacher to the \emph{thin teacher} by following the principle of \emph{hidden consistency} that will be introduced later. (ii) \textbf{Uniform Layer Selection} that selects layers of the thin teacher with the uniform strategy, and uses these parameters to initialize the student. (iii) \textbf{Transformer Distillation} that performs the last-layer distillation introduced before.

\begin{figure*}[t]
\includegraphics[width=\textwidth]{prune_process_v7.pdf}
\caption{Overview of width-wise extraction that includes the extraction of FFN neurons in the intermediate layer of FFN, head neurons in each attention head of MHA and the hidden neurons.}
\centering
\label{fig:pruneprocess}
\end{figure*}

\subsection{Width-wise Extraction.} The width-wise extraction includes the extraction of FFN neurons, head neurons and hidden neurons, which is illustrated in Fig.~\ref{fig:pruneprocess}. Specifically, we refer to the activation units in the intermediate layer of FFN as FFN neurons, which are represented by the gray circles in the figure, and call the activation units in each attention head of MHA as head neurons, which are represented by the green and blue circles in the figure. The hidden neurons, represented by the light blue circles in the figure, refer to the input and output activations of MHA and FFN, and the output of embedding layer.

Because there are residual connections in the MHA and FFN modules, we argue that the extraction of hidden neurons should follow the principle of {\it hidden consistency} which means that the extracted hidden neurons of different modules/layers should have the same position indexes. For example, as illustrated in Fig.~\ref{fig:pruneprocess}, the hidden neurons with some position indexes of 2 and 3 in different modules are extracted to meet the principle of hidden consistency. If the hidden consistency is not followed, the residual connection relationship will be destroyed, and the knowledge of reused teacher parameters will become chaotic. Experimental results in Table~\ref{tab:main} confirm that if the hidden consistency principle is violated, reusing parameters of the teacher model will even bring negative effects to the student model. In addition, we call the group of hidden neurons whose position index is $i$ in all sub-modules as the $i$-th hidden dimension. Extracting hidden neurons in our ETD strategy can be seen as selecting $d'_h$ hidden dimensions from $d_h$ ones, where $d'_h$ and $d_h$ refer to the hidden size of the student and the teacher respectively.

Following this principle, we propose two different approaches to extract the teacher's parameters. (i) \textbf{ETD-Rand}: We randomly extract neurons and assign the corresponding weight parameters to student model. (ii) \textbf{ETD-Impt}: We introduce the score-based pruning method in~\cite{sixteen} to extract the relatively important weights from teacher model. Specifically, in ETD-Impt, we calculate the importance of a neuron by the impact on the pre-training loss $\mathcal{L}$ if we remove it. Formally, we denote the $i$-th FFN neuron's value in the intermediate layer of FFN as $\bm{f_i}$. Using the first-order Taylor expansion, it’s importance $I_{f_i}$ can be calculated as:

\begin{equation}
\begin{split}
    I_{f_i} &= \left| \mathcal{L} - \mathcal{L}_{\bm{f_i}=0} \right|, \\
    &= \left| \mathcal{L} - (\mathcal{L}+ \frac{\partial \mathcal{L}}{\partial \bm{f_i}}(0-\bm{f_i}) + R_1(\bm{f_i}))\right|,\\
     &\approx \left|\frac{\partial \mathcal{L}}{\partial \bm{f_i}} \bm{f_i} \right|.
\end{split}
\end{equation}

The item $R_1(\bm{f_i})$ is the remainder of the first-order Taylor expansion so we can ignore it. Similarly, we denote the $i$-th head neuron’s value in a head as $\bm{h_i}$, whose importance $I_{h_i}$ can be calculated by:

\begin{equation}
    I_{h_i} = \left|\frac{\partial \mathcal{L}}{\partial \bm{h_i}}\bm{h_i} \right|.
\end{equation}

As shown in Fig.~\ref{fig:pruneprocess}, we extract the neurons in each head while keeping the number of attention heads. Because attention-based distillation has the constraint that the head number of student should be the same as the teacher. Then, we denote the $i$-th input hidden neuron of MHA module in the $l$-th Transformer layer as $n^l_i$ and the $j$-th hidden dimension as $d_j$. Their importance are denoted as $I_{n^l_i}$ and $I_{d_j}$ respectively, which can be calculated by:
\begin{align}
        \label{algorithn:prune_hidden}
        I_{n^l_i} &= \left|\frac{\partial \mathcal{L}}{\partial \bm{n^l_i}}\bm{n^l_i} \right|,    \\
        \label{algorithn:sum_hidden}
        I_{d_j} &= \sum_{l=0}^{N} I_{n^l_j},
\end{align}
where $\bm{n^l_i}$ means the value of the hidden neuron, $N$ means the number of the teacher's transformer layers. Our preliminary experiments found that summing the neurons’ importance of all Transformer layers in the same hidden dimension is a better choice to obtain $I_{d_i}$. Finally, after calculating the importance scores of all types of neurons, we extract neurons with the highest scores.

According to~\cite{sixteen}, we choose the first 10 thousand sentences in training data
to calculate the importance of neurons.\footnote{The FLOPs (floating-point operations) of calculating the importance score by 10 thousand sentences in training data is 2.5e15, we ignore this computational cost in Table~\ref{tab:main} because of the large difference of magnitude.}

\subsection{Uniform Layer Selection.} After width-wise extraction, we get a thin teacher model, which has the same width with the student. Assuming that the thin teacher and the student has $N$ and $M$ Transformer layers respectively, we choose $M$ layers from the thin teacher model to initialize the student model, and use the set $S$ to represent the chosen layers. It’s necessary to include 0 in $S$, because the index 0 refers to the embedding layer. Formally, the strategy of layer selection we apply is uniform-strategy ($S$ = $\left\{0,\lfloor\frac{N}{M}\rfloor,\lfloor\frac{2N}{M}\rfloor,...,N \right\}$). The other two typical strategies include top-strategy ($S$ = $\left\{0, N\!\!-\!\!M\!\!+\!\!1, N\!\! -\!\!M\!\!+\!\!2, ...,N\right\}$) and bottom strategy ($S$ = $\left\{0, 1, 2, ..., M\right\}$). Our experiments show that the uniform-strategy achieves better results than the other strategies.

\subsection{Transformer Distillation.} We empirically found that the last-layer distillation strategy performs better than the uniform one. Thus, we adopt the last-layer distillation strategy in ETD. As for the performance gap between these two strategies, we can see the experimental results shown in Table~\ref{tab:main} as ``ETD-TinyBERT" (uniform) and ``ETD-Impt" (last-layer).
\section{Experiments}
\subsection{Experimental Setup}

\noindent {\bf Datasets.} We use English Wikipedia and Toronto Book Corpus~\cite{book} for the distillation data. The maximum sequence length is 512. For evaluation, we use tasks from GLUE benchmark~\cite{wang2019glue} and SQuADv1.1~\cite{2016squad}. We report F1 for SQuADv1.1, Matthews correlation (Mcc) for CoLA and accuracy (Acc) for other tasks.

\noindent {\bf Implementation Details.} We use BERT-Base as our teacher. BERT-Base is a 12-layer Transformers with 768 hidden size. Our ETD method is generic and can be applied to students with varied architectures,to prove it, we instantiate two student models, which are 6-layer and 4-layer model with 384 hidden size respectively. They have different widths with the teacher.

We perform the last-layer distillation based on the released TinyBERT code\footnote{ https://github.com/huawei-noah/Pretrained-Language-Model/tree/master/TinyBERT} and follow the training setting: the batch size and the peak learning rate is set to 128 and 1e-4 respectively, the warm up proportion is set to 10\%. The training epochs of task-agnostic distillation is set to 5 by following~\cite{dynabert, layermapping}.\footnote{The training time that we distill a 6-layer student with 384 hidden size for 5 epochs is (80 hours * 8 V100 cards).} For GLUE fine-tuning, we set batch size to 32, choose the learning rate from \{5e-6, 1e-5, 2e-5, 3e-5\} and epochs from \{4, 5, 10\}. For SQuADv1.1 fine-tuning, we set batch size to 16, the learning rate to 3e-5 and the number of training epochs to 4.

\begin{table*}[t]
\centering
\resizebox{\textwidth}{!}{
\begin{tabular}{l|cc|cccccccc|c}
\hline
Strategy&\#Params&\#FLOPs&SST-2&MNLI&MRPC&CoLA&QNLI&QQP&STS-B&SQuADv1.1&Avg. \\
~ & ~ & (Training) & (Acc) & (Acc) & (Acc) & (Mcc) & (Acc) & (Acc) & (Acc) & (F1/EM) & ~ \\
\hline
BERT-Base(Teacher) &109M &6.43e19 &93.6 &84.7 &87.9 &59.6 &91.6 &91.4 &89.6  &88.4/81.0 &85.85 \\
DistilBERT\ddag\cite{distilbert} &66M & - & 91.3 & 82.2 &87.5 &51.3 &89.2 &88.5 & 86.9 &85.8/77.7 & 82.84 \\
\hline
\multicolumn{7}{l}{\it Main Results: Distill the 6-Layer Student Model with 384 Hidden Size}\\
\hline
Rand-Init       &22M &1.55e19 &90.9 &81.6 &86.0 &36.5&89.0&89.3&88.3&85.0/76.6 &80.83   \\
\textbf{ETD-Rand}   &22M &1.55e19 &91.6 &81.7 &\textbf{87.2} &37.7&89.2&89.2&88.1&85.2/76.8 &81.26   \\
\quad \textbf{43\% trained} &22M &6.64e18 &91.0 &81.6 &87.0 &35.9 &89.3 &\textbf{89.4} &88.2 &84.7/76.0 &80.88 \\
\textbf{ETD-Impt}   &22M &1.55e19 &91.3 &\textbf{82.0} &87.1 &\textbf{40.1}&\textbf{89.8}&89.2 & \textbf{88.5}&\textbf{85.6/77.4}& \textbf{81.70} \\
\quad \textbf{28\% trained} &22M &4.43e18 &\textbf{91.9} &81.4 &87.1 &37.2 &89.5 &89.2 &88.4  &84.7/76.1 &81.17 \\
ETD-wo-hidn-con & 22M &1.55e19 &90.6 &81.0 &86.0 &36.1 &89.1 &88.9 &88.2 & 84.9/76.1 & 80.60 \\
\hline
\multicolumn{7}{l}{\it Generality: Combine with Existing BERT Distillation Methods}\\
\hline
Ta-TinyBERT\dag\cite{tinybert}  &22M &1.55e19 &90.6 &80.9 &86.3 &34.9 &87.9 &88.6 &87.2  &83.9/75.2 &80.03 \\
\textbf{ETD-Ta-TinyBERT}  &22M &1.55e19 &91.1 &81.1 &86.0 &38.1 &88.3 &88.6 &87.6  &84.3/75.8&80.64\\
MiniLM\dag\cite{minilm}  &22M &1.55e19 &90.3 &81.4 &84.2 &35.7 &88.8 &89.1 &88.5 &84.5/75.7 &80.30\\
\textbf{ETD-MiniLM}  &22M &1.55e19 &91.1 &82.3 &86.5 &37.5 &89.4 &89.2 &88.4 &84.7/76.1 &81.14\\
\hline
\multicolumn{7}{l}{\it Generality: Distill the 4-Layer Student Model with 384 Hidden Size}\\
\hline 4L-Rand &19M&1.47e19&89.2 &79.2 &85.0 &25.4 & 87.9 &\textbf{88.5} &87.8 & 81.7/72.2 &78.07 \\
\textbf{4L-ETD-Impt} &19M&1.47e19&\textbf{89.6} &\textbf{79.4} &84.5 &\textbf{31.1} & \textbf{88.1} &88.3 &\textbf{88.0} & \textbf{81.9}/\textbf{72.4} &\textbf{78.86} \\
\quad \textbf{35\% trained}  &19M&5.25e18&89.4 &78.9 &\textbf{85.2} &27.4&87.5 &88.2 &87.8 &81.1/71.5 &78.18 \\
\hline
\end{tabular}
}
\caption{Comparison between the performance of student models with different initialization strategies. The fine-tuning results are all averaged over 3 runs on the dev set. \ddag\, denotes that the results are taken from \cite{mobilebert} and \dag\, means the models trained using the released code or the re-implemented code. We report the average score (Avg.) excluding the EM metric of SQuAD.}
\label{tab:main}
\end{table*}

\noindent {\bf Baselines.} As shown in Table~\ref{tab:main}, ``BERT-Base (Teacher)" is the teacher model of our distillation process. Additionally, we also report the results of DistilBERT~\cite{distilbert} as a baseline although it’s not a fair comparison because of the different model parameter amounts (DistilBERT is a 6-layer model whose width remains the same as the teacher model). ``Rand-Init" is our proposed baseline method which initializes the student randomly with the last-layer distillation strategy. The only difference with ETD is that Rand-Init does not reuse the teacher’s parameters. ``ETD-wo-hidn-con" is short for ``ETD-without-hidden-consistency", which denotes it violates the principle of hidden consistency. Specifically, we calculate the importance of the hidden neurons in each layer as shown in Eq.~(\ref{algorithn:prune_hidden}) and extract the most important ones inside each layer directly. Note that it does not guarantee that the position indexes of the extracted hidden neurons in each layer are the same, so it violates the principle. Moreover, in order to verify the generality of our ETD method, we apply ETD-Impt to different typical distillation methods (e.g., TinyBERT and MiniLM). ``TinyBERT" refers to TinyBERT which only does task-agnostic distillation in our paper, and ``ETD-TinyBERT/MiniLM" means that we use the ETD-Impt strategy to initialize the student model. We reproduced MiniLM ourselves because MiniLM doesn't share the model of this architecture setting.\footnote{The results differ from the original paper, possibly because we train with different hyper-parameters.} For a fair comparison, we train them with the same hyper-parameters.

\subsection{Main Results}

\begin{figure*}[t]
\centering
\includegraphics[width=0.9\textwidth]{gd_flops_loss_mnli_squad_v2.pdf}
\caption{The distillation loss curves and the score curves of MNLI and SQuADv1.1.}
\centering
\label{fig:gd_flops_loss}
\end{figure*}

The experimental results in and Table~\ref{tab:main} show that: 1) ETD-Rand/Impt strategies achieve comparable results to Rand-Init, while using only 43\% and 28\% computation cost (FLOPs)\footnote{https://github.com/google-research/electra/blob/master/flops computation.py} respectively. Compared with the baselines of TinyBERT and MiniLM, ETD-Impt can achieve the similar results with even less than 28\% computation cost. 2) ETD-Impt is better than ETD-Rand which indicates that the fine-grained extraction is needed to achieve better results. 3) Reusing the teacher’s parameters is beneficial for most tasks, especially for the tasks of CoLA and MRPC. 4) ETD is generic and can be combined with other distillation methods, such as TinyBERT or MiniLM, to improve their performance. Moreover, ETD works for the 4-layer architecture and achieves comparable results to Rand-Init with only 35\% of the training time. 5) We also find that ETD-Impt can perform comparably to DistilBERT in 7 out 8 tasks and is 3× smaller in size. 6) We demonstrate the effectiveness and necessity of following hidden consistency in Table~\ref{tab:main} and Fig.~\ref{fig:gd_flops_loss}. ``ETD-wo-hidn-con" poses a negative effect on the performance and performs even worse than the ``Rand-Init".

From the distillation loss shown in Fig.~\ref{fig:gd_flops_loss}, we can get the same conclusion that ETD method benefits the distillation efficiency. Moreover, the score curves of the two representative and complex tasks (MNLI and SQuADv1.1) are shown in Fig.~\ref{fig:gd_flops_loss}. We find that lower loss during distillation largely indicates better performance at downstream tasks.
\section{Ablation Study}
In this section, we conduct ablation experiments to analyze the effects of different proposed techniques on the performance, such as with-wise extraction and the strategy of layer selection. Note that all student models are 6-layer with 384 hidden size in this section.

\subsection{The Effect of Each Module}
To understand each type of neuron’s contribution, we also conduct the ablation experiments about the effect of FFN/head/hidden neurons in ETD-Impt. As shown in Table~\ref{tab:ablation}, ``ETD-Impt-rev-ffn" means that we extract the FFN neurons with the lowest importance scores instead of the highest ones. It can test the contribution of FFN neurons and the effectiveness of the importance score based method. It’s the same with ``ETD-Impt-rev-head" and ``ETD-Impt-rev-hidden". Specially, ``ETD-Impt-rev-all" means we reverse the ranks of all types of neurons. The results in Table~\ref{tab:ablation} shown that: 1) Extracting the least important neurons hurts the performance, which means that all three types of neurons have positive effects on the distillation process and the model’s performance in GLUE tasks and SQuAD. It also proves the effectiveness of our extracting method for these three types of neurons. Besides, we can know that among these three types of neurons, FFN neuron’s contribution is the smallest because the gap of distillation loss during convergence between ``ETD-Impt" and ``ETD-Impt-rev-ffn" is the closest. 2) The ETD method is always better than ``Rand-Init", even if it reuses the least important parameters as shown in the comparison between ``ETD-Imptrev-all" and ``Rand-Init", which further proves the effectiveness and necessity of reusing the teacher’s parameters.

\begin{table*}[t]
\centering
\resizebox{0.6\textwidth}{!}{
\begin{tabular}{l|c|cc|c}
\hline
Setting&Distillation Loss&Score&SQuADv1.1&Avg. \\
~ & (Convergence) & ~ & (F1/EM) & ~ \\
\hline
Rand-Init &1.231 & 80.24 & 85.0/76.6 & 80.83 \\
\hline \textbf{ETD-Impt} &\textbf{1.209}&\textbf{81.14}&\textbf{85.6/77.4}&\textbf{81.70}       \\
\hline ETD-Impt-rev-ffn &1.216&80.74&85.1/75.5 &81.29 \\
ETD-Impt-rev-head& 1.218&80.50&85.4/76.7&81.12  \\
ETD-Impt-rev-hidden& 1.219&80.66&85.0/76.6&81.20   \\
ETD-Impt-rev-all&1.227&80.38&84.9/76.3&80.95 \\
\hline
ETD-Impt-top& 1.215&80.55&85.3/76.8&81.14   \\
ETD-Impt-bottom&1.217&80.64&85.2/76.7&81.21 \\
\hline
\end{tabular}}
\caption{Ablation experiments to analyze the effects of different proposed techniques on the performance. The results are validated on the dev set and averaged over 3 runs. Distillation Loss refers to the final loss reached by the model at the end of training. Score refers to the average score of the GLUE tasks. We report the average score (Avg.) excluding the EM metric of SQuAD.}
\label{tab:ablation}
\end{table*}

\subsection{The Strategy of Layer Selection}
We compare the proposed uniform-strategy with the other two strategies: top-strategy (ETD-Impt-top) and bottom-strategy (ETD-Impt-bottom). According to the comparison results as shown in Table~\ref{tab:ablation}, we can find that: 1) The difference between bottom-strategy and top-strategy is not obvious. 2) The uniform strategy is the best performing strategy among them with obvious advantages. The results of the experiment indicate that uniform layer selection is a good strategy because it ensures that the student model can obtain the knowledge from the bottom to the top of the BERT model.
\section{Related Work}
\noindent {\bf Knowledge Distillation for BERT.} Knowledge distillation~\cite{kd} for BERT compression can be categorized into task-agnostic BERT distillation~\cite{distilbert,minilm, mobilebert} and task-specific BERT distillation~\cite{bertpkd, dynabert,tang2019lstm}. Task-agnostic distillation derives a small general BERT with pre-trained BERT as the teacher and large-scale unsupervised corpus as the training dataset. Task-specific BERT distillation learns a small fine-tuned BERT for a specific task with fine-tuned BERT as the teacher and task dataset as the input. TinyBERT~\cite{tinybert} proposes a two-stage distillation framework, including task-agnostic/specific distillation. Our work focuses on the task-agnostic distillation and aims to improve its efficiency.

\noindent {\bf Pruning.} Model extraction is essentially a pruning method. Han~\cite{prune} proposes an iterative pruning method based on weight magnitude. Similarly, the lottery tickets~\cite{ticket} use an iterative method to find the winning lottery ticket. However, these works both require specialized hardware optimizations for acceleration due to the irregular weight sparsity. There are some works of pruning-based Transformer compression~\cite{poorman, sixteen, l0, layerdrop, zhang2020know, dewynter2020optimal, gordon2020compressing}. In this paper, we propose the ETD strategy that adopts pruning techniques to extract FFN/head/hidden neurons of the teacher for finding a better initialization for the student.

\begin{comment}
Recently, Han et al. (2015) propose a
simple compression pipeline, achieving 40 times
reduction in model size without hurting accuracy.
Unfortunately, these techniques induce irregular
weight sparsity, which precludes highly optimized
computation routines. Thus, others explore pruning entire filters (Li et al., 2016; Liu et al., 2017),
with some even targeting device-centric metrics,
such as floating-point operations (Tang et al.,
2018) and latency (Chen et al., 2018). Still other
studies examine quantizing neural networks (Wu
et al., 2018); in the extreme, Courbariaux et al.
(2016) propose binarized networks with both binary weights and binary activations.

However, the high demand for computing resources in training such models hinders their application in practice. In order to alleviate this resource hunger in large-scale model training, we propose a Patient Knowledge Distillation approach to compress an original large model (teacher) into ....
\end{comment}

\section{Conclusion and Future Work}
In this paper, we proposed a flexible strategy of reusing teacher’s parameters named ETD for the efficient task-agnostic distillation. The experimental results show that our method significantly improves the efficiency of task-agnostic distillation and the fine-grained method ETD-Impt performs better than the random extraction ETD-Rand. In the future, we will optimize the strategy of layer selection and study more advanced pruning methods for ETD, such as L0 penalty~\cite{l0} and group lasso~\cite{grouplasso2006}.
\section{Example Appendix}

\end{document}