Asymmetry-plasma/main.tex

%% График от E^2 и E^4 в зависимости от концентрации (меняем eplion)
%%  FDTD Lumerical ?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This is the LaTeX ARTICLE template for RSC journals
%Copyright The Royal Society of Chemistry 2014
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\documentclass[twoside,twocolumn,9pt]{article}
\usepackage{extsizes}
\usepackage[super,sort&compress,comma]{natbib} 
\usepackage[version=3]{mhchem}
\usepackage[left=1.5cm, right=1.5cm, top=1.785cm, bottom=2.0cm]{geometry}
\usepackage{balance}
\usepackage{widetext}
\usepackage{times,mathptmx}
\usepackage{sectsty}
\usepackage{graphicx} 
\usepackage{lastpage}
\usepackage[format=plain,justification=raggedright,singlelinecheck=false,font={stretch=1.125,small,sf},labelfont=bf,labelsep=space]{caption}
\usepackage{float}
\usepackage{fancyhdr}
\usepackage{fnpos}
\usepackage{tabularx}
\usepackage[english]{babel}
\usepackage{array}
\usepackage{droidsans}
\usepackage{charter}
\usepackage[T1]{fontenc}
\usepackage[usenames,dvipsnames]{xcolor}
\usepackage{setspace}
\usepackage[compact]{titlesec}
%%% Please don't disable any packages in the preamble, as this may
%%% cause the template to display incorrectly.%%%
\usepackage{amsmath}
\newcommand{\red}[1]{{\color{red} #1}}

\usepackage{epstopdf}%This line makes .eps figures into .pdf - please
                     %comment out if not required.

\definecolor{cream}{RGB}{222,217,201}

\begin{document}

\pagestyle{fancy}
\thispagestyle{plain}
\fancypagestyle{plain}{

%%%HEADER%%%
\fancyhead[C]{\includegraphics[width=18.5cm]{head_foot/header_bar}}
\fancyhead[L]{\hspace{0cm}\vspace{1.5cm}\includegraphics[height=30pt]{head_foot/journal_name}}
\fancyhead[R]{\hspace{0cm}\vspace{1.7cm}\includegraphics[height=55pt]{head_foot/RSC_LOGO_CMYK}}
\renewcommand{\headrulewidth}{0pt}
}
%%%END OF HEADER%%%

%%%PAGE SETUP - Please do not change any commands within this
%%% section%%%
\makeFNbottom
\makeatletter
\renewcommand\LARGE{\@setfontsize\LARGE{15pt}{17}}
\renewcommand\Large{\@setfontsize\Large{12pt}{14}}
\renewcommand\large{\@setfontsize\large{10pt}{12}}
\renewcommand\footnotesize{\@setfontsize\footnotesize{7pt}{10}}
\makeatother

\renewcommand{\thefootnote}{\fnsymbol{footnote}}
\renewcommand\footnoterule{\vspace*{1pt}% 
\color{cream}\hrule width 3.5in height 0.4pt \color{black}\vspace*{5pt}} 
\setcounter{secnumdepth}{5}

\makeatletter 


\renewcommand\@biblabel[1]{#1}            
\renewcommand\@makefntext[1]% 
{\noindent\makebox[0pt][r]{\@thefnmark\,}#1}
\makeatother 
\renewcommand{\figurename}{\small{Fig.}~}
\sectionfont{\sffamily\Large}
\subsectionfont{\normalsize}
\subsubsectionfont{\bf}
\setstretch{1.125} %In particular, please do not alter this line.
\setlength{\skip\footins}{0.8cm}
\setlength{\footnotesep}{0.25cm}
\setlength{\jot}{10pt}
\titlespacing*{\section}{0pt}{4pt}{4pt}
\titlespacing*{\subsection}{0pt}{15pt}{1pt}
%%%END OF PAGE SETUP%%%

%%%FOOTER%%%
\fancyfoot{}
\fancyfoot[LO,RE]{\vspace{-7.1pt}\includegraphics[height=9pt]{head_foot/LF}}
\fancyfoot[CO]{\vspace{-7.1pt}\hspace{13.2cm}\includegraphics{head_foot/RF}}
\fancyfoot[CE]{\vspace{-7.2pt}\hspace{-14.2cm}\includegraphics{head_foot/RF}}
\fancyfoot[RO]{\footnotesize{\sffamily{1--\pageref{LastPage} ~\textbar  \hspace{2pt}\thepage}}}
\fancyfoot[LE]{\footnotesize{\sffamily{\thepage~\textbar\hspace{3.45cm} 1--\pageref{LastPage}}}}
\fancyhead{}
\renewcommand{\headrulewidth}{0pt} 
\renewcommand{\footrulewidth}{0pt}
\setlength{\arrayrulewidth}{1pt}
\setlength{\columnsep}{6.5mm}
\setlength\bibsep{1pt}
%%%END OF FOOTER%%%

%%%FIGURE SETUP - please do not change any commands within this
%%% section%%%
\makeatletter 
\newlength{\figrulesep} 
\setlength{\figrulesep}{0.5\textfloatsep} 

\newcommand{\topfigrule}{\vspace*{-1pt}% 
\noindent{\color{cream}\rule[-\figrulesep]{\columnwidth}{1.5pt}} }

\newcommand{\botfigrule}{\vspace*{-2pt}% 
\noindent{\color{cream}\rule[\figrulesep]{\columnwidth}{1.5pt}} }

\newcommand{\dblfigrule}{\vspace*{-1pt}% 
\noindent{\color{cream}\rule[-\figrulesep]{\textwidth}{1.5pt}} }

\makeatother
%%%END OF FIGURE SETUP%%%

%%%TITLE, AUTHORS AND ABSTRACT%%%
\twocolumn[
  \begin{@twocolumnfalse}
\vspace{3cm}
\sffamily
\begin{tabular}{m{4.5cm} p{13.5cm} }
  \includegraphics{head_foot/DOI} & \noindent\LARGE{\textbf{Photogenerated Electron-Hole Plasma Induced Symmetry Breaking in Spherical Silicon Nanoparticle}} \\%Article title goes here instead of the text "This is the title"
  \vspace{0.3cm} & \vspace{0.3cm} \\
                                  & \noindent\large{Anton Rudenko,$^{\ast}$\textit{$^{a}$}  Konstantin Ladutenko,\textit{$^{b}$}  Sergey Makarov\textit{$^{b}$} and Tatiana E. Itina\textit{$^{a}$$^{b}$}
 
                                    \textit{$^{a}$~Laboratoire Hubert Curien, UMR CNRS 5516, University of Lyon/UJM, 42000, Saint-Etienne, France }
                                    \textit{$^{b}$~ITMO University, Kronverksiy pr. 49, St. Petersburg, Russia}
 
                                    } \\%Author names go here instead of "Full name", etc.
 
 
  \includegraphics{head_foot/dates}
& \noindent\normalsize
  {The concept of nonlinear all-dielectric nanophotonics based on high
  refractive index (e.g., silicon) nanoparticles supporting magnetic
  optical response has recently emerged as a powerful tool for ultrafast
  all-optical modulation at nanoscale. A strong modulation can be
  achieved via photo-generation of dense electron-hole plasma in the
  regime of simultaneous excitation of electric and magnetic optical
  resonances, resulting in an effective transient reconfiguration of
  nanoparticle scattering properties. However, only homogeneous plasma
  generation was previously considered in the photo-excited
  nanoparticle, remaining unexplored any effects related to the
  plasma-induced optical inhomogeneities. Here we examine these
  effects by using 3D numerical modeling of coupled electrodynamic and
  material ionization equations. Based on the simulation results, we
  observed a deeply subwavelength plasma-induced nanopatterning of
  spherical silicon nanoparticles. In particular, we revealed strong
  symmetry breaking in the initially symmetrical nanoparticle, which
  arises during ultrafast photoexcitation near the magnetic dipole 
  resonance. The proposed ultrafast breaking of the nanoparticle
  symmetry paves the way to the novel opportunities for nonlinear
  optical nanodevices.} 

\end{tabular}



 \end{@twocolumnfalse} \vspace{0.6cm}

]%%%END OF TITLE, AUTHORS AND ABSTRACT%%%

%%%FONT SETUP - please do not change any commands within this section
\renewcommand*\rmdefault{bch}\normalfont\upshape \rmfamily
\section*{} \vspace{-1cm}


%%%FOOTNOTES%%%

\footnotetext{\textit{$^{a}$~Univ Lyon, UJM-St-Etienne, CNRS UMR 5516,
F-42000, Saint-Etienne, France}} \footnotetext{\textit{$^{b}$~ITMO
University, Kronverksiy pr. 49, St. Petersburg, Russia}}


% Please use \dag to cite the ESI in the main text of the article.

% If you article does not have ESI please remove the the \dag symbol
% from the title and the footnotetext below.
% \footnotetext{\dag~Electronic Supplementary Information (ESI)
% available: [details of any supplementary information available
% should be included here]. See DOI:10.1039/b000000x/} %additional
% addresses can be cited as above using the lower-case letters, c, d,
% e... If all authors are from the same address, no letter is required

% \footnotetext{\ddag~Additional footnotes to the title and authors can
% be included \emph{e.g.}\ `Present address:' or `These authors
% contributed equally to this work' as above using the symbols: \ddag,
% \textsection, and \P. Please place the appropriate symbol next to the
% author's name and include a \texttt{\textbackslash footnotetext} entry
% in the the correct place in the list.}


%%%END OF FOOTNOTES%%%

%%%MAIN TEXT%%%%
  
\section{Introduction}
 
All-dielectric nonlinear nanophotonics based on high refractive index
dielectric has become prospective paradigm in modern optics, owing to
recent advances in harmonics generation~\cite{shcherbakov2014enhanced,
  yang2015nonlinear, makarov2016self, shorokhov2016multifold,
  makarov2017efficient} and ultrafast all-optical
modulation~\cite{iyer2015reconfigurable, makarov2015tuning,
  shcherbakov2015ultrafast, yang2015nonlinear, baranov2016nonlinear,
  baranov2016tuning, shcherbakov2017ultrafast}. In fact,
all-dielectric nanoantennas and metasurfaces possess much smaller
parasitic Joule losses at high intensities as compared with their
plasmonic counterparts, whereas their nonlinear properties are
comparable. More importantly, the unique properties of the nonlinear
all-dielectric nanodevices are due to existing of both electric and
magnetic optical resonances in visible and near IR
ranges~\cite{kuznetsov2016optically}. For instance, even slight
variation of dielectric permittivity around optical resonances leads
to significant changes of optical properties (transmittance or
reflectance) of all-dielectric nanoantennas~\cite{makarov2015tuning,
  baranov2016nonlinear, baranov2016tuning} and
metasurfaces~\cite{iyer2015reconfigurable, shcherbakov2015ultrafast,
  yang2015nonlinear, shcherbakov2017ultrafast, makarov2017light}.

In these works on all-dielectric nonlinear nanostructures, the
building blocks (nanoparticles) were considered as objects with
dielectric permittivity \textit{homogeneously} distributed over
nanoparticle (NP). Therefore, in order to manipulate the propagation
angle of the transmitted light it was proposed to use complicated
nanostructures with reduced symmetry~\cite{albella2015switchable,
  baranov2016tuning, shibanuma2016unidirectional}.

\begin{figure}[t] \centering
  \includegraphics[width=0.75\linewidth]{Concept.pdf}
  \caption{Schematic illustration of electron-hole plasma 2D and 1D
distributions in silicon nanoparticle around a magnetic resonance.}
  \label{fgr:concept}
\end{figure}

Recently, plasma explosion imaging technique has been used to
observe electron-hole plasma (EHP), produced by femtosecond lasers,
inside NPs~\cite{Hickstein2014}. Particularly, a strongly
localized EHP in the front side\footnote{The incident wave propagates
  in positive direction of $z$ axis. For the NP with
  geometric center located at $z=0$ front side corresponds to the
  volume $z>0$ and back side for $z<0$} of NaCl nanocrystals of
$R = 100$~nm was revealed. The forward ejection of ions in this case
was attributed to a nanolensing effect inside the NP and the
intensity enhancement as low as $10\%$ on the far side of the
NP. Much stronger enhancements can be achieved near electric
and magnetic dipole resonances excited in single semiconductor
NPs, such as silicon (Si), germanium (Ge) etc.

In this Letter, we show that ultra-short laser-based EHP
photo-excitation in a spherical semiconductor (e.g., silicon) NP leads
to a strongly inhomogeneous carrier distribution. To reveal and study
this effect, we perform a full-wave numerical simulation. We consider
an intense femtosecond (\textit{fs}) laser pulse to interact with a
silicon NP supporting Mie resonances and two-photon free carrier
generation. In particular, we couple finite-difference time-domain
(FDTD) method used to solve three-dimensional Maxwell equations with
kinetic equations describing nonlinear EHP generation.
Three-dimensional transient variation of the material dielectric
permittivity is calculated for NPs of several sizes. The obtained
results propose a novel strategy to create complicated non-symmetrical
nanostructures by using single photo-excited spherical silicon
NPs. Moreover, we show that a dense EHP can be generated at deeply
subwavelength scale ($< \lambda / 10$) supporting the formation of
small metalized parts inside the NP. In fact, such effects transform
an all-dielectric NP to a hybrid metall-dielectric one strongly
extending functionality of the ultrafast optical nanoantennas.


%Plan:
%\begin{itemize}

%\item Fig.3: Temporal evolution of EHP in NP with fixed diameter (at
%MD) at different intensities, in order to show possible regimes of
%plasma-patterning of NP volume. It would be nice, if we will show
%power patterns decencies on intensity for side probe pulse to show beam steering due to symmetry breaking.
%\item Fig.4: (a) Dependence on pulse duration is also interesting. We
%have to show at which duration the asymmetry factor is saturated. (b)
%2D map of asymmetry factor in false colors, where x-axis and y-axis correspond to intensity and NP diameter.
%\end{itemize} %Additionally, if you will manage to calculate
%evolution of scattering power pattern and show considerable effect of
% beam steering, we can try Nanoscale or LPR, because the novelty will
% be very high.


\section{Modeling details}

We focus attention on silicon because this material is promising for
the implementation of numerous nonlinear photonic devices. This
advantage is based on a broad range of optical nonlinearities, strong
two-photon absorption, as well as a possibility of the photo-induced
EHP excitation~\cite{leuthold2010nonlinear}. Furthermore, silicon
nanoantennas demonstrate a sufficiently high damage threshold due to
the large melting temperature ($\approx 1690$~K), whereas its nonlinear
optical properties have been extensively studied during last
decays~\cite{Van1987, Sokolowski2000, leuthold2010nonlinear}. High
silicon melting point typically preserves structures formed from this
material up to the EHP densities on the order of the critical value
$N_{cr} \approx 5\cdot{10}^{21}$~cm$^{-3}$ \cite{Korfiatis2007}. At
the critical density and above, silicon acquires metallic properties
($Re(\epsilon) < 0$) and contributes to the EHP reconfiguration during
ultrashort laser irradiation.

The process of three-dimensional photo-generation of the EHP in
silicon NPs has not been modeled before in time-domain. Therefore,
herein we propose a model considering ultrashort laser interactions
with a resonant silicon sphere, where the EHP is generated via one-
and two-photon absorption processes.  Importantly, we also consider
nonlinear feedback of the material by taking into account the
intraband light absorption on the generated free carriers. To simplify
our model, we neglect free carrier diffusion due to the considered
short time scales. In fact, the aim of the present work is to study
the EHP dynamics \textit{during} ultra-short (\textit{fs}) laser
interaction with the NP. The created electron-hole modifies both
laser-particle interaction and, hence, the following particle
evolution. However, the plasma then will recombine at picosecond time
scale.

\subsection{Light propagation}

Ultra-short laser interaction and light propagation inside the silicon
NP are modeled by solving the system of three-dimensional Maxwell's equations
written in the following way
\begin{align} \begin{cases} \label{Maxwell}$$
    \displaystyle{\frac{\partial{\vec{E}}}{\partial{t}}=\frac{\nabla\times\vec{H}}{\epsilon_0\epsilon}-\frac{1}{\epsilon_0\epsilon}(\vec{J}_p+\vec{J}_{Kerr})} \\
    \displaystyle{\frac{\partial{\vec{H}}}{\partial{t}}=-\frac{\nabla\times\vec{E}}{\mu_0}},
    $$ \end{cases}
\end{align} 
where $\vec{E}$ is the electric field, $\vec{H}$ is the magnetizing
field, $\epsilon_0$ is the free space permittivity, $\mu_0$ is the
permeability of free space, $\epsilon = n^2 = 3.681^2$ is the
permittivity of non-excited silicon at $800$~nm wavelength
\cite{Green1995}, $\vec{J}_p$ and $\vec{J}_{Kerr}$ are the nonlinear
currents, which include the contribution due to Kerr effect
$\vec{J}_{Kerr} =
\epsilon_0\epsilon_\infty\chi_3\frac{\partial\left(\left|\vec{E}\right|^2\vec{E}\right)}{\partial{t}}$,
where $\chi_3 =4\cdot{10}^{-20}$ m$^2$/V$^2$ for laser wavelength
$\lambda = 800$~nm \cite{Bristow2007}, and heating of the conduction
band, described by the differential equation derived from the Drude
model
\begin{equation} \label{Drude}
  \displaystyle{\frac{\partial{\vec{J_p}}}{\partial{t}} = -
    \nu_e\vec{J_p} + \frac{e^2N_e(t)}{m_e^*}\vec{E}},
\end{equation}
where $e$ is the elementary charge, $m_e^* = 0.18m_e$ is the reduced
electron-hole mass \cite{Sokolowski2000}, $N_e(t)$ is the
time-dependent free carrier density and $\nu_e = 10^{15}$~s$^{-1}$ is
the electron collision frequency \cite{Sokolowski2000}. Silicon
NP is surrounded by vacuum, where the light propagation is
calculated by Maxwell's equations with $\vec{J} = 0$ and
$\epsilon = 1$. The system of Maxwell's equations coupled with
electron density equation is solved by the finite-difference numerical
method \cite{Rudenko2016}, based on the finite-difference time-domain
(FDTD) \cite{Yee1966} and auxiliary-differential methods for
dispersive media \cite{Taflove1995}. At the edges of the grid, we
apply the absorbing boundary conditions related to convolutional
perfectly matched layers (CPML) to avoid nonphysical reflections
\cite{Roden2000}. The initial electric field is introduced as a
Gaussian slightly focused beam as follows
\begin{align}
\begin{aligned}
\label{Gaussian}
{E_x}(t, r, z) = \frac{w_0}{w(z)}{exp}\left(i\omega{t} - \frac{r^2}{{w(z)}^2} - ikz - ik\frac{r^2}{2R(z)} + i\varsigma(z)\right)\\
\times\;{exp}\left(-\frac{4\ln{2}(t-t_0)^2}{\theta^2}\right),
\end{aligned}
\end{align}
where $\theta \approx 80$~\textit{fs} is the temporal pulse width at the half maximum (FWHM),
$t_0$ is a time delay, $w_0 = 3{\mu}m$ is the waist beam,
$w(z) = {w_0}\sqrt{1+(\frac{z}{z_R})^2}$ is the Gaussian's beam spot
size, $\omega = 2{\pi}c/{\lambda}$ is the angular frequency,
$\lambda = 800$~nm is the laser wavelength in air, $c$ is the speed of
light, ${z_R} = \frac{\pi{{w_0}^2}n_0}{\lambda}$ is the Rayleigh
length, $r = \sqrt{x^2 + y^2}$ is the radial distance from the beam's
waist, $R_z = z\left(1+(\frac{z_R}{z})^2\right)$ is the radius of
curvature of the wavelength comprising the beam, and
$\varsigma(z) = {\arctan}(\frac{z}{z_R}) $ is the Gouy phase shift.

\subsection{Material ionization}

To account for the material ionization that is induced by a
sufficiently intense laser field inside the particle, we couple
Maxwell's equations with the kinetic equation for the electron-hole
plasma as described below.

The time-dependent conduction-band carrier density evolution is
described by a rate equation that was proposed by van Driel
\cite{Van1987}. This equation takes into account such processes as
photoionization, avalanche ionization and Auger recombination, and is
written as
\begin{equation} \label{Dens}
  \displaystyle{\frac{\partial{N_e}}{\partial t} =
    \frac{N_a-N_e}{N_a}\left(\frac{\sigma_1I}{\hbar\omega} +
      \frac{\sigma_2I^2}{2\hbar\omega}\right) + \alpha{I}N	_e -
    \frac{C\cdot{N_e}^3}{C\tau_{rec}N_e^2+1},} \end{equation} where
$I=\frac{n}{2}\sqrt{\frac{\epsilon_0}{\mu_0}}\left|\vec{E}\right|^2$
is the intensity, $\sigma_1 = 1.021\cdot{10}^3$~cm$^{-1}$ and
$\sigma_2 = 0.1\cdot{10}^{-7}$~cm/W are the one-photon and two-photon
interband cross-sections \cite{Choi2002, Bristow2007, Derrien2013},
$N_a = 5\cdot{10}^{22}$~cm$^{-3}$ is the saturation particle density
\cite{Derrien2013}, $C = 3.8\cdot{10}^{-31}$~cm$^6$/s is the Auger
recombination rate \cite{Van1987}, $\tau_{rec} = 6\cdot{10}^{-12}$~s is
the minimum Auger recombination time \cite{Yoffa1980}, and
$\alpha = 21.2$~cm$^2$/J is the avalanche ionization coefficient
\cite{Pronko1998} at the wavelength $800$~nm in air. As we have noted,
free carrier diffusion is neglected during and shortly after the laser
excitation \cite{Van1987, Sokolowski2000}. In particular, from the
Einstein formula $D = k_B T_e \tau/m^* \approx (1$--$\,2)\cdot{10}^{-3}$ m$^2$/s
($k_B$ is the Boltzmann constant, $T_e$ is the electron temperature,
$\tau=1$~\textit{fs} is the collision time, $m^* = 0.18 m_e$ is the effective
mass), where $T_e \approx 2*{10}^4$ K for $N_e$ close to $N_{cr}$ \cite{Ramer2014}. It
means that during the pulse duration ($\approx 80$~\textit{fs}) the diffusion
length will be around 5$\,$--10~nm for $N_e$ close to $N_{cr}$.

\begin{figure}[ht!] 
\centering
\includegraphics[width=0.495\textwidth]{mie-fdtd-3}
\caption{\label{mie-fdtd} (a, b) First four Lorentz-Mie coefficients
  ($a_1$, $a_2$, $b_1$, $b_2$) and factors of asymmetry $G_I$,
  $G_{I^2}$ according to Mie theory at fixed wavelength $800$~nm. (c,
  d) Squared intensity distribution calculated by Mie theory and (e,
  f) EHP distribution for low free carrier densities
  $N_e \approx 10^{20}$~cm$^{-3}$ by Maxwell's equations
  (\ref{Maxwell}, \ref{Drude}) coupled with EHP density equation
  (\ref{Dens}). (c-f) Incident light propagates from the left to the
  right along $Z$ axis, electric field polarization $\vec{E}$ is along
  $X$ axis.}
\end{figure}

The changes of the real and imaginary parts of the permittivity
associated with the time-dependent free carrier response
\cite{Sokolowski2000} can be derived from equations (\ref{Maxwell},
\ref{Drude}) and are written as follows
\begin{align} \begin{cases}  \label{Index} $$
	\displaystyle{Re(\epsilon^\ast) = \epsilon -\frac{{e^2}N_e}{\epsilon_0m_e^*(\omega^2+{\nu_e}^2)}} \\
	\displaystyle{Im(\epsilon^\ast) = \frac{{e^2}N_e\nu_e}{\epsilon_0m_e^*\omega(\omega^2+{\nu_e}^2)}.}
 $$ \end{cases} \end{align}

\subsection{Mie calculations}
A steady-state interaction of a plain electromagnetic wave with a
spherical particle has a well-known analytical solution described by a
Mie theory~\cite{Bohren1983}. It is only valid in the absence of
nonlinear optical response, thus we can compare it against
above-mentioned FDTD-EHP model only for small plasma densities, where
we can neglect EHP impact to the refractive index. Non-stationary
nature of a femtosecond pulse increase the complexity of the
analysis. A detailed discussion on the relation between Mie theory and
FDTD-EHP model will be provided in the next section.

We used Scattnlay program to evaluate calculations of Mie coefficients
and near-field distribution~\cite{Ladutenko2017}.  This program is
available online at GitHub~\cite{Scattnlay-web} under open source
license.

\section{Results and discussion}

\begin{figure*}[p]
 \centering
 \includegraphics[width=145mm]{time-evolution-I-no-NP.pdf}
 \caption{\label{time-evolution} Temporal EHP (a, c, e) and asymmetry
   factor $G_{N_e}$ (b, d, f) evolution for different Si nanoparticle
   radii of (a, b) $R = 75$~nm, (c, d) $R = 100$~nm, and (e, f)
   $R = 115$~nm. Pulse duration $80$~\textit{fs}
   (FWHM). \red{\textbf{TODO:} on the plot it looks more than 100 fs
     for FWHM!!! Anton? } Wavelength $800$~nm in air. (b, d, f)
   Different stages of EHP evolution shown in Fig.~\ref{plasma-grid}
   are indicated. The temporal evolution of the incident Gaussian beam
   intensity is also shown. Peak laser fluence is fixed to be
   $0.125$~J/cm$^2$.}
\vspace*{\floatsep}
 \centering
 \includegraphics[width=150mm]{plasma-grid.pdf}
 \caption{\label{plasma-grid} EHP density snapshots inside Si nanoparticle of
   radii $R = 75$~nm (a-d), $R = 100$~nm (e-h) and $R = 115$~nm (i-l)
   taken at different times and conditions of excitation (stages
   $1-4$: (1) first optical cycle, (2) extremum at few optical cycles,
   (3) Mie theory, (4) nonlinear effects). $\Delta{Re(\epsilon)}$
   indicates the real part change of the dielectric function defined
   by Equation (\ref{Index}). Pulse duration $80$~\textit{fs}
   (FWHM). Wavelength $800$~nm in air. Peak laser fluence is fixed to
   be $0.125$~J/cm$^2$.}
 \end{figure*}

%\subsection{Effect of the irradiation intensity on EHP generation}

 We start with a pure electromagnetic problem without EHP
 generation. We plot in Fig.~\ref{mie-fdtd}(a) Mie coefficients of a
 Si NP as a function of its size for a fixed laser wavelength
 $\lambda = 800$~nm.  For the NP sizes under consideration most of
 contribution to the electromagnetic response originates from electric
 and magnetic dipole (ED and MD), while for sizes near
 $R \approx 150$~nm a magnetic quadrupole (MQ) response turns to be
 the main one. The superposition of multipoles defines the
 distribution of electric field inside of the NP.  We introduce $G_I$
 factor of asymmetry, corresponding to difference between the volume
 integral of intensity in the front side of the NP to that in the back
 side normalized to their sum:
 $G_I = (I^{front}-I^{back})/(I^{front}+I^{back})$, where
 $I^{front}=\int_{(z>0)}|E|^2d{\mathrm{v}}$ and
 $I^{back}=\int_{(z<0)} |E|^2d{\mathrm{v}}$ are expressed via
 amplitude of the electric field $|E|$. The factor $G_{I^2}$ was
 determined in a similar way by using volume integrals of squared
 intensity to predict EHP asymmetry due to two-photon absorption.
 Fig.~\ref{mie-fdtd}(b) shows $G$ factors as a function of the NP
 size. For the NPs of sizes below the first MD resonance,
 the intensity is enhanced in the front side as in
 Fig.~\ref{mie-fdtd}(c) and $G_I > 0$. The behavior changes near the
 size resonance value, corresponding to $R \approx 105$~nm. In
 contrast, for larger sizes, the intensity is enhanced in the back
 side of the NP as demonstrated in Fig.~\ref{mie-fdtd}(d). In fact,
 very similar EHP distributions can be obtained by applying Maxwell's
 equations coupled with the rate equation for relatively weak
 excitation with EHP concentration of $N_e \approx 10^{20}$~cm$^{-3}$,
 see Fig.~\ref{mie-fdtd}(e,f). The optical properties do not change
 considerably due to the excitation according to
 (\ref{Index}). Therefore, the excitation processes follow the
 intensity distribution. However, such coincidence was achieved under
 quasi-stationary conditions, after the electric field made enough
 oscillations inside the Si NP. Further on we present transient
 analysis, which reveals much more details.

 To achieve a quantative description for evolution of the EHP
 distribution during the \textit{fs} pulse, we introduced another
 asymmetry factor
 $G_{N_e} = (N_e^{front}-N_e^{back})/(N_e^{front}+N_e^{back})$
 indicating the relationship between the average EHP densities in the
 front and in the back halfs of the NP. This way, $G_{N_e} = 0$
 corresponds to the quasi-homogeneous case and the assumption of the
 NP homogeneous EHP distribution can be made to investigate the
 optical response of the excited Si NP. When $G_{N_e}$ significantly
 differs from $0$, this assumption, however, could not be
 justified. In what follows, we discuss the results of the numerical
 modeling (see Fig.~\ref{time-evolution}) of the temporal evolution of
 EHP densities and the asymmetry factor $G_{N_e}$. It reveals the EHP
 evolution stages during pulse duration. Typical change of the
 permittivity corresponding to each stage is shown in
 Fig.~\ref{plasma-grid}.  For better visual representation of time
 scale of the whole incident pulse and its single optical cycle we put a
 squared electric field profile on all plots in
 Fig.~\ref{time-evolution} in gray color as a backgroud image (note
 linear time scale on the left column and logarithmic scale on the
 right one).

 To describe all the stages of light non-linear interaction with Si
 NP, we present the calculation results obtained by using Maxwell's
 equations coupled with electron kinetics equations for different
 radii for resonant and non-resonant conditions. In this case, the
 geometry of the EHP distribution can strongly deviate from the
 intensity distribution given by Mie theory. Two main reasons cause
 the deviation: (i) non-stationarity of interaction between
 electromagnetic pulse and NP (ii) nonlinear effects, taking place due
 to transient optical changes in Si. The non-stationary intensity
 deposition during \textit{fs} pulse results in different time delays
 for exciting electric and magnetic resonances inside Si NP because of
 different quality factors $Q$ of the resonances.

 In particular, MD resonance (\textit{b1}) has $Q \approx 8$, whereas
 electric one (\textit{a1}) has $Q \approx 4$. The larger particle
 supporting MQ resonance (\textit{b2}) demonstrates $ Q \approx
 40$. As soon as the electromagnetic wave period at $\lambda = 800$~nm
 is 2.7~\textit{fs}, one needs about 10~\textit{fs} to pump the ED,
 20~\textit{fs} for the MD, and about 100~\textit{fs} for the MQ.
 According to these considerations, after few optical cycles taking
 place on a 10~\textit{fs} scale it results in the excitation of the
 low-\textit{Q} ED resonance, which dominates MD and MQ independently
 on the exact size of NPs. Moreover, during the first optical cycle
 there is no multipole modes structure inside of NP, which results
 into a very similar field distribution for all size of NP under
 consideration as shown in Fig.~\ref{plasma-grid}(a,e,i) . We address
 to this phenomena as \textit{'Stage~1'}. This stage demonstrates the
 initial penetration of electromagnetic field into the NP during the
 first optical cycle.
 
 \textit{'Stage~2'} corresponds to further electric field oscillations
 ($t \approx 5$--$15$) leading to the formation of ED field pattern in
 the center of the NP as it can be seen in
 Fig.~\ref{plasma-grid}(f,j). We would like to stress the
 unstationery nature of field pattern at this stage. The energy
 balance between extinction and pumping is not set, moreover, there is
 a simultaneous growth of the incident pulse apmlitude. This leads to
 a superposition of ED field pattern with the one from the Stage 1,
 resulting into the presence for the maximum of the EHP distribution
 in the front side of the Si NP. This effect dominates for the
 smallest NP with $R=75$~nm in Fig.~\ref{plasma-grid}(b), where ED
 mode is tuned far away from the resonance (see Fig.~\ref{mie-fdtd}(c)
 for field supression inside of NP predicted by Mie theory). At this
 stage, the density of EHP ($N_e < 10^{20}$~cm$^2$) is still not high
 enough to significantly affect the optical properties of the NP.

 When the number of optical cycles is big enough ($t>20$~\textit{fs})
 both ED and MD modes can be exited to the level necessary to achieve
 the stationary intensity pattern corresponding to the Mie-based
 intensity distribution at the \textit{'Stage~3'} (see
 Fig.~\ref{plasma-grid}). The EHP density for the most volume of NP is
 still relatively small to affect the EHP evolution or for diffusion,
 but is already high enough to change the local optical
 properties. Below the MD resonance $R \approx 100$~nm, the EHP is
 mostly localized in the front side of the NP as shown in
 Fig.~\ref{plasma-grid}(c). The highest stationary asymmetry factor
 $G_{N_e} \approx 0.5$--$0.6$ is achieved in this case. At the MD
 resonance conditions, the EHP distribution has a toroidal shape and
 is much closer to the homogeneous distribution. In contrast, above
 the MD resonant size for $R = 115$~nm, and the $G_{N_e} < 0$ due to
 the fact that EHP is dominantly localized in the back side of the NP.

 Once again, due to presence of continous pumping the Stage~3 is
 superposed with Stage~1 field pattern, resulting in the EHP localized
 in the front side. This can be seen when comparing result from the
 Mie theory in Fig.~\ref{mie-fdtd}(d) and result of full 3D simulation
 in Fig.~\ref{mie-fdtd}(f). Note that pumping of NP significantly
 changes during a single optical cycle, this leads to a large
 variation of asymmetry factor $G_{N_e}$ at first stage. This
 variation stedialy decrease as it goes to Stage~3.

 The explain this we need to consider time evolution of mean EHP
 densities $N_e$ in the front and back halfs of NP presented in
 Fig.~\ref{time-evolution}(a,c,e). As soon as recombination and
 diffusion processes are negligible at \textit{fs} time scale, both
 $N_e^{front}$ and $N_e^{back}$ curves experience monotonous behaviour
 with small pumping steps synced to the incident pulse. Front and back
 halfs are separated in space, wich obviously leads to the presence of
 time delay between puping steps in each curve caused with the same
 optical cycle of the incident wave. This delay causes a large value
 of the assymetry factor during first stage. However, as soon as mean
 EHP density increases the contribution of this pumping steps to
 resulting assymetry becomes smallar and the variation of $G_{N_e}$
 synced with the period of incident light decreases.

 %A bookmark by Kostya
 
 For the higher excitation conditions, the optical properties of
 silicon change significantly according to the equations
 (\ref{Index}). As a result, the non-resonant ED
 contributes to the forward shifting of EHP density
 maximum. Therefore, EHP is localized in the front part of the NP,
 influencing the asymmetry factor $G_{N_e}$ in
 Fig.~\ref{time-evolution}. Approximately at the pulse peak, the critical
 electron density $N_{cr} = 5\cdot{10}^{21}$~cm$^{-3}$ for silicon,
 which corresponds to the transition to quasi-metallic state
 $Re(\epsilon) \approx 0$ and to the electron plasma resonance, is
 overcome. Further irradiation leads to a decrease in the asymmetry
 parameter down to $G_{N_e} = 0$ for higher EHP densities, as one can
 observe in Fig.~\ref{plasma-grid}(d, h, l).

 As the EHP acquires quasi-metallic properties at stronger excitation
 $N_e > 5\cdot{10}^{21}$~cm$^{-3}$, the EHP distribution evolves
 inside NPs because of the photoionization and avalanche ionization
 induced transient optical response and the effect of newly formed
 EHP. This way, the distribution becomes more homogeneous and the
 effect is likely to be enhanced by electron diffusion inside Si
 NPs. We refer to these nonlinear phenomena as \textit{'Stage~4'}.
     
 It is worth noting that it is possible to achieve a formation of
 deeply subwavelength EHP regions due to high field localization. The
 smallest EHP localization and the larger asymmetry factor are
 achieved below the MD resonant conditions for $R < 100$~nm.
 Thus, the EHP distribution in Fig.~\ref{plasma-grid}(c) is optimal for
 symmetry breaking in Si NP, as it results in the larger asymmetry
 factor $G_{N_e}$ and higher electron densities $N_e$. We stress here
 that such regime could be still safe for NP due to the very small
 volume where such high EHP density is formed.

% \subsection{Effects of NP size and scattering efficiency
% factor on scattering directions}

% \begin{figure}[ht] \centering
% \includegraphics[width=90mm]{time-evolution.png}
% \caption{\label{time-evolution} a) Scattering efficiency factor $Q_{sca}$
% dependence on the radius $R$ of non-excited silicon NP
% calculated by Mie theory; b) Parameter of forward/backward scattering
% dependence on the radius $R$ calculated by Mie theory for non-excited
% silicon NP c) Optimization parameter $K$ dependence on the
% average electron density $n_e^{front}$ in the front half of the
% NP for indicated radii (1-7).}
% \end{figure}

% We have discussed the EHP kinetics for a silicon NP of a
% fixed radius $R \approx 105$~nm. In what follows, we investigate the
% influence of the NP size on the EHP patterns and temporal
% evolution during ultrashort laser irradiation. A brief analysis of
% the initial intensity distribution inside the NP given by
% the classical Mie theory for homogeneous spherical particles
% \cite{Mie1908} can be useful in this case. Fig. \ref{time-evolution}(a, b)
% shows the scattering efficiency and the asymmetry parameter for
% forward/backward scattering for non-excited silicon NPs of
% different radii calculated by Mie theory \cite{Mie1908}. Scattering
% efficiency dependence gives us the value of resonant sizes of
% NPs, where the initial electric fields are significantly
% enhanced and, therefore, we can expect that the following conditions
% will result in a stronger electron density gradients. Additionally,
% in the case of maximum forward or backward scattering, the initial
% intensity distribution has the maximum of asymmetry. One can note,
% that for $R \approx 100$~nm and $R \approx 150$~nm both criteria are
% fulfilled: the intensity is enhanced $5-10$ times due to
% near-resonance conditions and its distribution has a strong
% asymmetry.

% In what follows, we present the calculation results obtained by
% using Maxwell's equations coupled with electron kinetics for
% different extremum radii for resonant and non-resonant
% conditions. One can note, that the maximum asymmetry factor of EHP
% $G$ does not guarantee the optimal asymmetry of intensity
% distribution, as the size of generated plasma and the value of the
% electron density equally contribute to the change of the modified
% NP optical response. For example, it is easier to localize
% high electron densities inside smaller NPs, however, due
% to the negligible size of the generated EHP with respect to laser
% wavelength in media, the intensity distribution around the
% NP will not change considerably. Therefore, we propose to
% introduce the optimization factor
% $K = \frac{n_e(G-1)R^2}{n_{cr}G{R_0^2}}$, where $R_0 = 100$~nm,
% $n_{cr} = 5\cdot{10}^{21}~cm^{-3}$, and $G$ is asymmetry of EHP,
% defined previously. The calculation results for different radii of
% silicon NPs and electron densities are presented in
% Fig. \ref{time-evolution}(c). One can see, that the maximum value are achieved
% for the NPs, that satisfy both initial maximum forward
% scattering and not far from the first resonant condition. For larger
% NPs, lower values of EHP asymmetry factor are obtained, as
% the electron density evolves not only from the intensity patterns in
% the front side of the NP but also in the back side.

%TODO: 
%Need to discuss agreement/differences between Mie and FDTD+rate. Anton ?!

% To demonstrate the effect of symmetry breaking, we calculate the
% intensity distribution around the NP for double-pulse
% experiment. The first pulse of larger pulse energy and polarization
% along $Ox$ generates asymmetric EHP inside silicon NP,
% whereas the second pulse of lower pulse energy and polarization $Oz$
% interacts with EHP after the first pulse is gone. The minimum
% relaxation time of high electron density in silicon is
% $\tau_{rec} = 6\cdot{10}^{-12}$~s \cite{Yoffa1980}, therefore, the
% electron density will not have time to decrease significantly for
% subpicosecond pulse separations. In our simulations, we use
% $\delta{t} = 200\:f\!s$ pulse separation. The intensity
% distributions near the silicon NP of $R = 95$~nm,
% corresponding to maxima value of $K$ optimization factor, without
% plasma and with generated plasma are shown in Fig. \ref{fig4}. The
% intensity distribution is strongly asymmetric in the case of EHP
% presence. One can note, that the excited NP is out of
% quasi-resonant condition and the intensity enhancements in
% Fig. \ref{fig4}(c) are weaker than in Fig. \ref{fig4}(b). Therefore,
% the generated nanoplasma acts like a quasi-metallic nonconcentric
% nanoshell inside the NP, providing a symmetry reduction
% \cite{Wang2006}.

% \begin{figure}[ht] \centering
% \includegraphics[width=90mm]{fig4.png}
% \caption{\label{fig4} a) Electron plasma distribution inside Si
% NP $R \approx 95$~nm 50~\textit{fs} after the pulse peak;
% (Kostya: Is it scattering field intensity snapshot $XXX\:f\!s$ after
% the second pulse maxima passed the particle?) Intensity
% distributions around and inside the NP b) without plasma,
% c) with electron plasma inside.}
% \end{figure}

%\begin{figure} %\centering
% \includegraphics[height=0.7\linewidth]{Si-flow-R140-XYZ-Eabs}
% \caption{EHP distributions for nonres., MD, ED, and MQ NPs
% at moderate photoexcitation. The aim is to show different possible
% EHP patterns and how strong could be symmetry breaking.
% \label{fgr:example}
%\end{figure}

%\subsection{Asymmetry analysis: effects of pulse duration, intensity
% and size} It is important to optimize asymmetry by varying pulse
% duration, intensity and size.

\section{Conclusions} We have considered ultra-short and sufficiently
intense light interactions with a single semiconductor nanoparticle
under different irradiation conditions and for various particle
sizes. As a result of the presented self-consistent calculations, we
have obtained spatio-temporal EHP evolution inside the
NPs and investigated the asymmetry of EHP distributions. % for different laser intensities. % and temporal pulse widths.
%It has been demonstrated that the EHP generation strongly affects
%NP scattering and, in particular, changes the preferable
%scattering direction. 
Different pathways of EHP evolution from the front side to the back
side have been revealed, depending on the NP sizes, and the
origins of different behavior have been explained by the
non-stationarity of the energy deposition and different quality
resonant factors for exciting the electric and magnetic dipole
resonances, intensity distribution by Mie theory and newly
plasma-induced nonlinear effects. The effect of the strong broadband
electric dipole resonance on the EHP asymmetric distribution during
first optical cycles has been revealed for different size
parameters. The higher EHP asymmetry is established for NPs
of smaller sizes below the first magnetic dipole
resonance. Essentially different EHP evolution and lower asymmetry is
achieved for larger NPs due to the stationary intensity
enhancement in the back side of the NP. The EHP densities
above the critical value were shown to lead to the EHP distribution
homogenization.
% In particular, the scattering efficiency factor is used to define
% the optimum NP size for preferential forward or backward
% scattering. Furthermore, a parameter has been introduced to describe
% the scattering asymmetry as a ratio of the EHP density in the front
% side to that in the back side of the NP. This parameter
% can be then used for two-dimensional scattering mapping, which is
% particularly important in numerous photonics applications.
The EHP asymmetry opens a wide range of applications in NP
nanomashining/manipulation at nanoscale, in catalysis as well as numerous
nano-bio-applications.  The observed plasma-induced breaking symmetry
can be also useful for beam steering, or for the enhanced second
harmonics generation.

\section{Acknowledgments}
A. R. and T. E. I. gratefully acknowledge the CINES of CNRS for
computer support. S. V. M. is thankful to ITMO Fellowship
Program. This work was partially supported by Russian Foundation for
Basic Researches (grants 17-03-00621, 17-02-00538, 16-29-05317).




%%%END OF MAIN TEXT%%%

%The \balance command can be used to balance the columns on the final
%page if desired. It should be placed anywhere within the first column
%of the last page.

%\balance

%If notes are included in your references you can change the title
% from 'References' to 'Notes and references' using the following
% command:

% \renewcommand\refname{Notes and references}

%%%REFERENCES%%%

\bibliography{References} %You need to replace "rsc" on this line
                              %with the name of your .bib file
\bibliographystyle{rsc} %the RSC's .bst file

\end{document}