reflow

main.tex

@@ -241,26 +241,34 @@ 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}. On the other hand, plasma explosion imaging technique~\cite{Hickstein2014} revealed \textit{in situ} strongly asymmetrical electron-hole plasma (EHP) distribution in various dielectric NPs during their pumping by femtosecond laser pulses. Therefore, local permittivity in the strongly photoexcited NPs can be significantly inhomogeneous, and symmetry of nanoparticles can be reduced. 
+  baranov2016tuning, shibanuma2016unidirectional}. On the other hand,
+plasma explosion imaging technique~\cite{Hickstein2014} revealed
+\textit{in situ} strongly asymmetrical electron-hole plasma (EHP)
+distribution in various dielectric NPs during their pumping by
+femtosecond laser pulses. Therefore, local permittivity in the
+strongly photoexcited NPs can be significantly inhomogeneous, and
+symmetry of nanoparticles can be reduced.
 
 %The forward ejection of ions was attributed in this case to a nanolensing effect inside the NP and to 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 theoretically that ultra-fast photo-excitation in a spherical silicon NP leads
-to a strongly inhomogeneous EHP distribution, as schematically shown in Fig.~\ref{fgr:concept}. To reveal and analyze 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 EHP
-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.
+In this Letter, we show theoretically that ultra-fast photo-excitation
+in a spherical silicon NP leads to a strongly inhomogeneous EHP
+distribution, as schematically shown in Fig.~\ref{fgr:concept}. To
+reveal and analyze 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 EHP 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
-a dielectric NP to a hybrid metall-dielectric one strongly
-extending functionality of the ultrafast optical nanoantennas.
+small metalized parts inside the NP. In fact, such effects transform a
+dielectric NP to a hybrid metall-dielectric one strongly extending
+functionality of the ultrafast optical nanoantennas.
 
 
 %Plan:
@@ -297,24 +305,29 @@ 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 and temporal evolution of EHP in
-silicon NPs has not been modeled before. 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.
+The process of three-dimensional photo-generation and temporal
+evolution of EHP in silicon NPs has not been modeled
+before. 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}
 
-The incident wave propagates in positive direction of $z$ axis, the NP geometric center located at $z=0$ front side corresponds to the volume $z>0$ and back side for $z<0$, as shown in Fig.~\ref{fgr:concept}. 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
+The incident wave propagates in positive direction of $z$ axis, the NP
+geometric center located at $z=0$ front side corresponds to the volume
+$z>0$ and back side for $z<0$, as shown in
+Fig.~\ref{fgr:concept}. 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}},
@@ -358,10 +371,10 @@ Gaussian slightly focused beam as follows
 \times\;{exp}\left(-\frac{4\ln{2}(t-t_0)^2}{\theta^2}\right),
 \end{aligned}
 \end{align}
-where $\theta \approx 130$~\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,
+where $\theta \approx 130$~\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
@@ -384,7 +397,7 @@ 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{\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
@@ -392,32 +405,35 @@ $\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
+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 130$~\textit{fs}) the diffusion
-length will be around 10$\,$--15~nm for $N_e$ close to $N_{cr}$.
+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 130$~\textit{fs}) the diffusion length will be around
+10$\,$--15~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) First four Lorentz-Mie coefficients
   ($a_1$, $a_2$, $b_1$, $b_2$) and (b) factors of asymmetry $G_I$,
-  $G_{I^2}$ according to the Mie theory at fixed wavelength $800$~nm. (c,
-  d) Squared intensity distributions at different radii \textit{R} calculated by the Mie theory and (e,
-  f) EHP distribution for low densities
-  $N_e \approx 10^{20}$~cm$^{-3}$ by Maxwell's equations
-  (\ref{Maxwell}, \ref{Drude}) coupled with EHP density equation
-  (\ref{Dens}). \red{Please, write $G_I$ and $G_{I^2}$ instead of \textit{G} as a title to axis-Y in (b)} (c-f) Incident light propagates from the left to the
-  right along $Z$ axis, electric field polarization $\vec{E}$ is along
-  $X$ axis.}
+  $G_{I^2}$ according to the Mie theory at fixed wavelength
+  $800$~nm. (c, d) Squared intensity distributions at different radii
+  \textit{R} calculated by the Mie theory and (e, f) EHP distribution
+  for low densities $N_e \approx 10^{20}$~cm$^{-3}$ by Maxwell's
+  equations (\ref{Maxwell}, \ref{Drude}) coupled with EHP density
+  equation (\ref{Dens}). \red{Please, write $G_I$ and $G_{I^2}$
+    instead of \textit{G} as a title to axis-Y in (b)} (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
@@ -437,13 +453,13 @@ 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~\textit{fs} pulse increases the complexity of the
-analysis. A detailed discussion on the relation between the Mie theory and
-FDTD-EHP model will be provided in the next section.
+analysis. A detailed discussion on the relation between the Mie theory
+and FDTD-EHP model will be provided in the next section.
 
-We used the Scattnlay program to evaluate calculations of Lorentz-Mie coefficients ($a_i$, $b_i$)
-and near-field distribution~\cite{Ladutenko2017}. This program is
-available online at GitHub~\cite{Scattnlay-web} under open source
-license.
+We used the Scattnlay program to evaluate calculations of Lorentz-Mie
+coefficients ($a_i$, $b_i$) and near-field
+distribution~\cite{Ladutenko2017}. This program is available online at
+GitHub~\cite{Scattnlay-web} under open source license.
 
 \section{Results and discussion}
 
@@ -452,28 +468,34 @@ license.
 \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 $130$~\textit{fs} (FWHM).  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
-  Gaussian beam intensity is also shown. Peak laser intensity is fixed to
-   be 10$^{12}$~W/cm$^2$. For better visual representation of time scale at a single optical cycle we put a  squared electric field profile in all plots in  Fig.~\ref{time-evolution} in gray color as a background image (note  linear time scale on the left column and logarithmic scale on the  right one). \red{Please, remove italic style of font for 't, fs' on titles of X-axises. Also, all fonts should be the same: please don't mix Times and Arial.}}
+ \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 $130$~\textit{fs} (FWHM).  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
+   Gaussian beam intensity is also shown. Peak laser intensity is
+   fixed to be 10$^{12}$~W/cm$^2$. For better visual representation of
+   time scale at a single optical cycle we put a squared electric
+   field profile in all plots in Fig.~\ref{time-evolution} in gray
+   color as a background image (note linear time scale on the left
+   column and logarithmic scale on the right one). \red{Please, remove
+     italic style of font for 't, fs' on titles of X-axises. Also, all
+     fonts should be the same: please don't mix Times and Arial.}}
     \end{figure*}
     
 \begin{figure*}   
  \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) maximum during few optical cycles,
-   (3) quasi-stationary regime, (4) strongly nonlinear regime). $\Delta{Re(\epsilon)}$
-   indicates the real part change of the local permittivity defined
-   by Equation (\ref{Index}). Pulse duration 130~fs
-   (FWHM). Wavelength 800~nm in air. Peak laser intensity is fixed to
-   be 10$^{12}$~W/cm$^2$. }
+ \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) maximum
+   during few optical cycles, (3) quasi-stationary regime, (4)
+   strongly nonlinear regime). $\Delta{Re(\epsilon)}$ indicates the
+   real part change of the local permittivity defined by Equation
+   (\ref{Index}). Pulse duration 130~fs (FWHM). Wavelength 800~nm in
+   air. Peak laser intensity is fixed to be 10$^{12}$~W/cm$^2$. }
  \end{figure*}
 
 %\subsection{Effect of the irradiation intensity on EHP generation}
@@ -488,8 +510,8 @@ license.
  the main one. The superposition of multipoles defines the
  distribution of electric field inside the NP. We introduce $G_I$
  factor of asymmetry, corresponding to difference between the volume
- integral of squared electric field in the front side ($I^{front}$) of the NP to that in the back
- side ($I^{back}$) normalized to their sum:
+ integral of squared electric field in the front side ($I^{front}$) of
+ the NP to that in the back side ($I^{back}$) 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
@@ -497,127 +519,145 @@ license.
  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,
- rather 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.
+ 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, rather 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 quantitative 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 halves of the NP, defined as $N_e^{front}=\frac{2}{V}\int_{(z>0)} {N_e}d{\mathrm{v}}$ and
- $N_e^{back}=\frac{2}{V}\int_{(z<0)} {N_e}d{\mathrm{v}}$, where $V = \frac{4}{3}\pi{R}^3$ is the volume of the nanosphere. In this way, $G_{N_e} = 0$ corresponds to the symmetrical 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 of the temporal evolution of
- EHP densities and the asymmetry factors $G_{N_e}$ for different sizes of Si NP, as shown in Fig.~\ref{time-evolution}). Typical change of the
- permittivity corresponding to each stage is shown in
- Fig.~\ref{plasma-grid}. It reveals the EHP
- evolution stages during interaction of femtosecond laser pulse with the Si NPs. 
+ indicating the relationship between the average EHP densities in the
+ front and in the back halves of the NP, defined as
+ $N_e^{front}=\frac{2}{V}\int_{(z>0)} {N_e}d{\mathrm{v}}$ and
+ $N_e^{back}=\frac{2}{V}\int_{(z<0)} {N_e}d{\mathrm{v}}$, where
+ $V = \frac{4}{3}\pi{R}^3$ is the volume of the nanosphere. In this
+ way, $G_{N_e} = 0$ corresponds to the symmetrical 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 of the temporal evolution of EHP densities and the
+ asymmetry factors $G_{N_e}$ for different sizes of Si NP, as shown in
+ Fig.~\ref{time-evolution}). Typical change of the permittivity
+ corresponding to each stage is shown in Fig.~\ref{plasma-grid}. It
+ reveals the EHP evolution stages during interaction of femtosecond
+ laser pulse with the Si NPs.
  
  \subsection{Stages of transient Si nanoparticle photoexcitation}
 
-To describe all the stages of non-linear light interaction with Si
+ To describe all the stages of non-linear light 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 at peak intensity 10$^{12}$~W/cm$^2$ and $\lambda = 800~nm$. In this case, the
- geometry of the EHP distribution can strongly deviate from the
- intensity distribution given by the Mie theory. Two main reasons cause
- the deviation: (i) non-stationarity of interaction between
+ radii for resonant and non-resonant conditions at peak intensity
+ 10$^{12}$~W/cm$^2$ and $\lambda = 800~nm$. In this case, the geometry
+ of the EHP distribution can strongly deviate from the intensity
+ distribution given by the 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 $\approx 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.
+ 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
+ $\approx 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 NP, which results
- in a very similar field distribution for all sizes of NP under
+ there is no multipole modes structure inside NP, which results in a
+ very similar field distribution for all sizes 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. Resulting factors $G_{N_e}$ exhibit sharp increase at this stage (Fig.~\ref{time-evolution}(b,d,f), yielding strong and ultrafast symmetry breaking. 
+ first optical cycle. Resulting factors $G_{N_e}$ exhibit sharp
+ increase at this stage (Fig.~\ref{time-evolution}(b,d,f), yielding
+ strong and ultrafast symmetry breaking.
  
  \textit{'Stage~2'} corresponds to further electric field oscillations
- ($t \approx 5$--$15$~fs) leading to the formation of ED E-field pattern in
- the center of the Si NP as shown in
- Fig.~\ref{plasma-grid}(f,j). We stress the
- nonstationary nature of E-field pattern at this stage, whereas there is
- a simultaneous growth of the incident pulse amplitude. This leads to
- a superposition of ED near-field pattern with that from the Stage 1,
- resulting in EHP concentration 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 suppression inside NP predicted by the 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 (Figs.~\ref{time-evolution}(a,c,e)).
-
- When the number of optical cycles is large 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}(c,g,k)). The EHP density for the most volume of NP is
- still relatively small to affect the EHP evolution,
- but is already high enough to change the local optical
- properties, i.e. real part of permittivity. Below the MD resonance ($R = 75$~nm), the EHP is
- mostly localized in the front side of the NP as shown in
- Fig.~\ref{plasma-grid}(c). The highest quasi-stationary asymmetry factor
- $G_{N_e} \approx 0.5$--$0.6$ is achieved in this case (Fig~\ref{time-evolution}(b)). At the MD
- resonance conditions ($R = 100$~nm), 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 the $G_{N_e} < 0$ due to
- the fact that EHP is dominantly localized in the back side of the NP.
-
- In other words, due to a quasi-stationary pumping during the Stage~3 is
- superposed with the Stage~1 field pattern, resulting in an additional
- 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 steadily decrease as it goes to Stage~3, as shown in Fig.~\ref{time-evolution}(b,d,f).
+ ($t \approx 5$--$15$~fs) leading to the formation of ED E-field
+ pattern in the center of the Si NP as shown in
+ Fig.~\ref{plasma-grid}(f,j). We stress the nonstationary nature of
+ E-field pattern at this stage, whereas there is a simultaneous growth
+ of the incident pulse amplitude. This leads to a superposition of ED
+ near-field pattern with that from the Stage 1, resulting in EHP
+ concentration 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 suppression inside NP predicted by
+ the 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
+ (Figs.~\ref{time-evolution}(a,c,e)).
+
+ When the number of optical cycles is large 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}(c,g,k)). The EHP density for the most
+ volume of NP is still relatively small to affect the EHP evolution,
+ but is already high enough to change the local optical properties,
+ i.e. real part of permittivity. Below the MD resonance ($R = 75$~nm),
+ the EHP is mostly localized in the front side of the NP as shown in
+ Fig.~\ref{plasma-grid}(c). The highest quasi-stationary asymmetry
+ factor $G_{N_e} \approx 0.5$--$0.6$ is achieved in this case
+ (Fig~\ref{time-evolution}(b)). At the MD resonance conditions
+ ($R = 100$~nm), 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 the $G_{N_e} < 0$ due to the fact that
+ EHP is dominantly localized in the back side of the NP.
+
+ In other words, due to a quasi-stationary pumping during the Stage~3
+ is superposed with the Stage~1 field pattern, resulting in an
+ additional 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 steadily decrease as it goes to Stage~3,
+ as shown in Fig.~\ref{time-evolution}(b,d,f).
 
  To explain this effect, we consider the time evolution of average EHP
  densities $N_e$ in the front and back halves of NP presented in
  Fig.~\ref{time-evolution}(a,c,e). As soon as the recombination and
  diffusion processes are negligible at \textit{fs} time scale, both
  $N_e^{front}$ and $N_e^{back}$ curves experience monotonous behavior
- with small pumping steps synced to the incident pulse. The front and the back
- halves of NP are separated in space, which obviously leads to the presence of
- time delay between pumping steps in each curve caused by the same
- optical cycle of the incident wave. This delay causes a large asymmetry factor during first stage. However, as soon as average
- EHP density increases the relative contribution of this pumping steps to
- the resulting asymmetry becomes smaller. This way variations of asymmetry
- $G_{N_e}$ synced with the period of incident light decreases.
+ with small pumping steps synced to the incident pulse. The front and
+ the back halves of NP are separated in space, which obviously leads
+ to the presence of time delay between pumping steps in each curve
+ caused by the same optical cycle of the incident wave. This delay
+ causes a large asymmetry factor during first stage. However, as soon
+ as average EHP density increases the relative contribution of this
+ pumping steps to the resulting asymmetry becomes smaller. This way
+ variations of asymmetry $G_{N_e}$ synced with the period of incident
+ light decreases.
 
  Higher excitation conditions are followed by larger values of
  electric field amplitude, which lead to the appearance of high EHP
  densities causing a significant change in the optical properties of
- silicon according to the equations (\ref{Index}). From the Mie theory, the initial (at the end of Stage~3) spatial pattern of the
- optical properties is non-homogeneous. When non-homogeneity of
-the optical properties becomes strong enough it leads to the
- reconfiguration of the E-field inside NP, which in turn strongly affects further reconfiguration of the optical properties. We
- refer to these strong nonlinear phenomena as \textit{'Stage~4'}. In
+ silicon according to the equations (\ref{Index}). From the Mie
+ theory, the initial (at the end of Stage~3) spatial pattern of the
+ optical properties is non-homogeneous. When non-homogeneity of the
+ optical properties becomes strong enough it leads to the
+ reconfiguration of the E-field inside NP, which in turn strongly
+ affects further reconfiguration of the optical properties. We refer
+ to these strong nonlinear phenomena as \textit{'Stage~4'}. In
  general, the reconfiguration of the electric field is unavoidable as
  far as the result from the Mie theory comes with the assumption of
  homogeneous optical properties in a spherical NP.
@@ -625,14 +665,18 @@ the optical properties becomes strong enough it leads to the
  Thus, the evolution of EHP density during Stage~4 depends on the
  result of multipole modes superposition at the end of Stage~3 and is
  quite different as we change the size of NP. For $R=75$~nm and
- $R=100$~nm, we observe a front side asymmetry before Stage~4, however,
- the origin of it is quite different. The $R=75$~nm NP is out of
- resonance, moreover, Mie field pattern and the one, which comes from
- the Stage~1 are quite similar. As soon as EHP density becomes high enough
- to change optical properties, the NP is still out of resonance,
- however, the presence of EHP increases absorption in agreement with
- (\ref{Index}), decreases Q-factor, and destroys optical modes.
- %This effect effectively leads to a partial screening, and it  becomes harder for the incident wave to penetrate deeper into EHP. Finally,  this finishes spilling the NP`s volume with plasma reducing the asymmetry, see Fig.~\ref{plasma-grid}(d).
+ $R=100$~nm, we observe a front side asymmetry before Stage~4,
+ however, the origin of it is quite different. The $R=75$~nm NP is out
+ of resonance, moreover, Mie field pattern and the one, which comes
+ from the Stage~1 are quite similar. As soon as EHP density becomes
+ high enough to change optical properties, the NP is still out of
+ resonance, however, the presence of EHP increases absorption in
+ agreement with (\ref{Index}), decreases Q-factor, and destroys
+ optical modes.
+ % This effect effectively leads to a partial screening, and it
+ % becomes harder for the incident wave to penetrate deeper into
+ % EHP. Finally, this finishes spilling the NP`s volume with plasma
+ % reducing the asymmetry, see Fig.~\ref{plasma-grid}(d).
  
  For $R=100$~nm, the evolution during the final stage goes in a
  similar way, with a notable exception regarding MD resonance.  As
@@ -643,20 +687,20 @@ the optical properties becomes strong enough it leads to the
  mark.
  
  The last NP with $R=115$~nm shows the most complex behavior during
- the Stage~4. The superposition of Mie-like E-field pattern with that from
- Stage~1 results in the presence of two EHP spatial maxima, back and
- front shifted. They serve as starting seeds for the EHP formation,
- and an interplay between them forms a complex behavior of the asymmetry
- factor curve. Namely, the sign is changed from negative to positive
- and back during the last stage. This numerical result can hardly be
- explained in a simple qualitative manner, it is too complex to
- account all near-field interaction of incident light with two EHP
- regions inside a single NP. It is interesting to note, however, that in
- a similar way as it was for $R=100$~nm the increased absorption
- should destroy ED and MD resonances, which are responsible for the back-shifted
- EHP. As soon as this EHP region is quite visible on the last snapshot
- in Fig.~\ref{plasma-grid}(l), this means that EHP seeds are
- self-supporting.
+ the Stage~4. The superposition of Mie-like E-field pattern with that
+ from Stage~1 results in the presence of two EHP spatial maxima, back
+ and front shifted. They serve as starting seeds for the EHP
+ formation, and an interplay between them forms a complex behavior of
+ the asymmetry factor curve. Namely, the sign is changed from negative
+ to positive and back during the last stage. This numerical result can
+ hardly be explained in a simple qualitative manner, it is too complex
+ to account all near-field interaction of incident light with two EHP
+ regions inside a single NP. It is interesting to note, however, that
+ in a similar way as it was for $R=100$~nm the increased absorption
+ should destroy ED and MD resonances, which are responsible for the
+ back-shifted EHP. As soon as this EHP region is quite visible on the
+ last snapshot in Fig.~\ref{plasma-grid}(l), this means that EHP seeds
+ are self-supporting.
  
  %A bookmark by Kostya
 
@@ -784,11 +828,25 @@ the optical properties becomes strong enough it leads to the
 % and size} It is important to optimize asymmetry by varying pulse
 % duration, intensity and size.
 
-\section{Conclusions} We have rigorously modeled and studied ultra-fast and intense light interaction with a single silicon nanoparticle of various sizes for the first time to our best knowledge. As a result of the presented self-consistent nonlinear calculations, we have obtained spatio-temporal EHP evolution inside the
-NPs and investigated the asymmetry of the EHP distributions. 
-We have revealed EHP strong asymmetric distribution during first optical cycles for different sizes. The highest average EHP asymmetry has been observed for NPs of smaller sizes below the first magnetic dipole resonance, when EHP is concentrated in the front side mostly during the laser pulse absorption. Essentially different EHP evolution and lower asymmetry has been achieved for larger NPs due to the  intensity enhancement in the back side of the NP. The EHP densities
-above the critical value have been shown to lead to homogenization of the EHP distribution.
-The observed plasma-induced breaking symmetry can be useful for creation of nonsymmetrical nanophotonic designes, e.g. for beam steering or enhanced second harmonics generation. Also, the asymmetric EHP opens a wide range of applications in NP nanomashining at deeply subwavelength scale. 
+ \section{Conclusions} We have rigorously modeled and studied
+ ultra-fast and intense light interaction with a single silicon
+ nanoparticle of various sizes for the first time to our best
+ knowledge. As a result of the presented self-consistent nonlinear
+ calculations, we have obtained spatio-temporal EHP evolution inside
+ the NPs and investigated the asymmetry of the EHP distributions.  We
+ have revealed EHP strong asymmetric distribution during first optical
+ cycles for different sizes. The highest average EHP asymmetry has
+ been observed for NPs of smaller sizes below the first magnetic
+ dipole resonance, when EHP is concentrated in the front side mostly
+ during the laser pulse absorption. Essentially different EHP
+ evolution and lower asymmetry has been achieved for larger NPs due to
+ the intensity enhancement in the back side of the NP. The EHP
+ densities above the critical value have been shown to lead to
+ homogenization of the EHP distribution.  The observed plasma-induced
+ breaking symmetry can be useful for creation of nonsymmetrical
+ nanophotonic designes, e.g. for beam steering or enhanced second
+ harmonics generation. Also, the asymmetric EHP opens a wide range of
+ applications in NP nanomashining at deeply subwavelength scale.
 
 \section{Acknowledgments}
 A. R. and T. E. I. gratefully acknowledge the CINES of CNRS for