# Relaxation to steady states of a binary liquid mixture around an optically
heated colloid

Takeaki Araki  Department of Physics, Kyoto University, Kyoto 606-8502, Japan.
Juan Ruben Gomez-Solano  Instituto de Física, Universidad Nacional Autónoma de
Mexico, Ciudad de Mexico, C.P. 04510, Mexico  Anna Maciołek
<EMAIL_ADDRESS>Institute of Physical Chemistry, Polish Academy of
Sciences, Kasprzaka 44/52, PL-01-224 Warsaw, Poland  Max-Planck-Institut für
Intelligente Systeme Stuttgart, Heisenbergstr. 3, D-70569 Stuttgart, Germany

(September 3, 2024)

###### Abstract

We study the relaxation dynamics of a binary liquid mixture near a light-
absorbing Janus particle after switching on and off illumination using
experiments and theoretical models. The dynamics is controlled by the
temperature gradient formed around the heated particle. Our results show that
the relaxation is asymmetric: the approach to a nonequilibrium steady state is
much slower than the return to thermal equilibrium. Approaching a
nonequilibrium steady state after a sudden temperature change is a two-step
process that overshoots the response of spatial variance of the concentration
field. The initial growth of concentration fluctuations after switching on
illumination follows a power law in agreement with the hydrodynamic and purely
diffusive model. The energy out-flow from the system after switching off
illumination is well described by a stretched exponential function of time
with characteristic time proportional to the ratio of the energy stored in the
steady state to the total energy flux in this state.

##  I Introduction

Relaxation processes are fundamental phenomena in soft and condensed matter
physics. The most research attention has been paid to relaxation processes
close to equilibrium, which are now well understood Onsager_1931a ;
Onsager_1931b ; kubo1957statistical ; dattagupta2012relaxation ;
metzler1999anomalous ; shiraishi2019information . This is not the case for
relaxation near nonequilibrium steady states, although various aspects of
transient nonequilibrium dynamics have been addressed in recent studies
PhysRevLett.125.110602 ; PhysRevLett.107.010601 ; Gomes_2017 ;
PhysRevE.99.032132 ; PhysRevLett.103.040601 ; PhysRevE.79.060104 ;
Gieseler_2015 ; PhysRevX.10.011066 ; PhysRevLett.120.180604 . One of the basic
questions here is what determines the timescales for a given system to reach a
nonequilibrium steady state and to relax back to equilibrium. This problem is
of particular relevance for light-absorbing particles that heat up under laser
illumination of an appropriate wavelength baffou2013thermo . In response to an
emerging temperature gradient, the surrounding medium restructures giving rise
to a variety of curious phenomena, which can be of great practical use. For
example, irradiated gold nanoparticles enhance locally a lipid membrane
permeability Rossi:2017 , which is used in various photothermal therapies, hot
nanoparticles can trigger cell fusion Oddershede:2018 , optically heated
colloids may self-propel Volpe-et:2011 ; Buttinioni-et:2012 ;
PhysRevLett.116.138301 ; 10.3389/fphy.2020.570842 ; vutukuri2020light ;
auschra2021thermotaxis  or, if trapped by a laser beam, they become
microscopic engines girot2016motion ; Volpe-et:2018 ; schmidt2021non .
Switching on and off illumination is the way to steer these processes, which
is of particular relevance when the fluid medium surrounding the particle
phase-separates in response to a temperature variation onuki2002phase . Then,
it is of immediate interest to know how quickly the medium near the particle
changes its structure after switching on illumination, compared with the speed
of return to the equilibrium state after switching it off and how this can be
controlled. Understanding relaxations of binary liquid mixtures in presence of
temperature gradients is also of greatest relevance for applications, such as
in optical fluid pumps doi:10.1063/5.0015247 , phase separation in
microfluidic cavities PhysRevFluids.6.024001 , and phoretic motion of colloids
in liquid mixtures PhysRevResearch.2.033177 .

Here, we address these issues in a system, which over the last decade has
become a paradigm for studying various aspects of active matter, i.e., a Janus
silica particle half-coated by gold and suspended in a binary liquid mixture
RevModPhys.88.045006 . Initially, the binary solvent at its critical
composition is prepared in the mixed state at temperature $T_{0}=28^{\circ}$C
below its lower critical point $T_{c}$. After applying green laser
illumination onto the sample- see Fig.1(a), the temperature non-isotropically
increases around the particle surface. This is due to the absorption peak of
gold around $\lambda=532$ nm and the poor absorption of silica and of the
surrounding fluid at that wavelength. With sufficient laser intensity, the
induced temperature quench exceeds locally $T_{c}$. This results in local
demixing of a binary solvent around the colloid, which is responsible for its
active motion PhysRevE.97.042603 ; C8SM01258J ; PhysRevLett.115.188304 ;
gomez2017tuning ; C9SM00509A .

Figure 1: (a) Schematic representation of the side-view ($x$-$z$ plane) of the
sample cell composed of two glass plates at separation $h\approx 2R$, where
$R$ is the particle radius, which is suspended in a PnP-water mixture at its
critical composition below $T_{c}$. At such confinement the particle is
immobilized. Green laser illumination of intensity $I$ is applied onto the
sample in the $z$-direction from above. (b) Sketch of the main particle axis
along which the spatial variance of the order parameters is computed up to
distance $L$ from the particle surface. The dark area on the right side of the
colloid corresponds to the golden cap.

In the nonequilibrium steady state one observes a stable droplet covering the
Janus particle asymmetrically; the droplet is much more pronounced near the
hot golden cap (see Fig.1(b)). Our experimental data and results of numerical
calculations within two theoretical models reveal that the reorganization of
the concentration field around the Janus colloid until the steady state is
reached is remarkably slow, whereas returning to the initial homogeneous
distribution after switching off illumination is relatively fast. This
asymmetry has not been noticed so far. The relaxation to a steady state and
back to equilibrium in our system has curious features due to the presence of
a temperature gradient. Among them is the power law of the early relaxation
after switching on illumination (ON process). The power-law exponent $\approx
2$, as predicted from both purely diffusive and hydrodynamic models, is in
excellent agreement with the experiment. This agreement supports the existence
of a generic mechanism of the initial process connected to the formation of
surface layers D0SM00964D . On the other hand, the relaxation to equilibrium
after switching off the illumination (OFF process) is stretched exponential
with a characteristic time decreasing as the illumination power increases,
which we attribute to the temperature gradient previously formed by the laser-
induced heating of the particle’s golden cap.

##  II Experimental results

In our experiments, we use a binary liquid mixture of propylene glycol
n-propyl ether (PnP) and water, the lower critical point (LCP) of which is
$T_{c}=31.9^{\circ}$C at its critical composition (0.4 PnP mass fraction)
Bauduin2004 . In order to characterize the relaxation of such a binary liquid
around the heated colloid, we study the time evolution of the image intensity
profile $-\frac{1}{2}<\phi<\frac{1}{2}$ measured along the main particle axis,
where $\phi$ is related to the concentration profile D0SM00964D . The value
$\phi=0$ corresponds to the fully mixed fluid, whereas $\phi<0$ and $\phi>0$
represent locally water-rich and PnP-rich regions, respectively. Spatio-
temporal variations of $\phi$ originate from the difference between the
refractive index of water and PnP within the temperature range of the
experiments ($30\pm 5^{\circ}$C), $n=1.331\pm 0.001$ and $n=1.410\pm 0.002$,
respectively. In the mixed phase at $T_{0}<T_{c}$ and without laser
illumination, the refractive index of the binary liquid in thermal equilibrium
is homogeneous, thereby leading to a uniform light intensity pattern recorded
by the camera around the Janus colloid. On the other hand, if the power of the
incident laser is sufficiently high such that the golden cap reaches the
temperature $T>T_{c}$, spatial variations of the refractive index due to the
local phase separation of the fluid give rise to deviations of the light
intensity recorded by the camera from the uniform intensity under thermal
equilibrium conditions, which allows us to implement shadowgraph visualization
ASSENHEIMER1994373 ; Mauger .

At a given time $t$, we compute the spatial variance of $\phi$ along
$R<r<R+L$, defined as

$\langle\Delta\phi(t)^{2}\rangle=\frac{1}{L}\int_{R}^{R+L}dr\,\left[\phi(r,t)-\langle\phi(t)\rangle\right]^{2},$
(1)

where

$\langle\phi(t)\rangle=\frac{1}{L}\int_{R}^{R+L}dr\,\phi(r,t),$ (2)

is the spatial average of the order parameter $\phi(r,t)$ at time $t$ along
the main particle axis, as shown in Fig.1(b). For all the calculations
presented in the following, $L=25\,\mu$m, which is approximately two times the
particle radius. Note that when the concentration profile has reached a steady
state, $\langle\Delta\phi(t)^{2}\rangle$ remains constant over time, whereas a
time dependence of $\langle\Delta\phi(t)^{2}\rangle$ reveals the transient
behavior of the concentration field of the fluid during ON and OFF processes.
On the other hand, the minimum value of $\langle\Delta\phi(t)^{2}\rangle$ is
set by the intensity fluctuations on the colormap for a fully mixed binary
liquid mixture, whose standard deviation is $\delta\phi_{0}\approx 0.03$.

Figure 2: Time evolution of the spatial variance of the order parameter
$\phi(r,t)$, computed along the main particle axis $R<r<R+L$ (with
$L=25\,\mu$m), for different values of the temperature increase $\Delta T$
induced by optical heating. Panels (a), (b), (c) and (d) correspond to $\Delta
T=7.8,9.1,10.4$ and $11.7^{\circ}$C, respectively. At time $t=-10$ s the laser
illumination is switched on, thus leading the formation a nonequilibrium
steady concentration field around the particle, while at $t=0$ it is switched
off. Insets: images of the steady state at $t=-1$ s. Figures 2(e) and 2(f)
show the initial growth of the spatial variance from the experiment ($\Delta
T=7.8{{}^{\circ}}$C) and from the numerical simulations ($h=4$) of the
hydrodynamic model, respectively. $\delta\phi_{0}^{2}=0.001$ and
$\delta\varphi_{0}^{2}=0.00303$

Figure 2 shows the time evolution of $\langle\Delta\phi(t)^{2}\rangle$ for
different values of the temperature increase $\Delta
T=T-T_{0}=7.8,9.1,10.4,11.7^{\circ}$C of the golden cap of a Janus particle
(radius $R=12\,\mu$m, strongly hydrophilic cap and weakly hydrophilic silica
hemisphere) above the bath temperature $T_{0}=28^{\circ}\mathrm{C}<T_{c}$ and
$T>T_{c}$ the temperature of the golden cap under green illumination. At time
$t=-10$ s, the illumination is switched on, which leads to the transient
coarsening process investigated in our previous paper D0SM00964D . Due to the
heating of the particle, the concentration profile changes with respect to the
initial uniform concentration in thermal equilibrium at temperature $T_{0}$,
for which $\langle\phi(t)\rangle=0$ and
$\langle\Delta\phi(t)^{2}\rangle=\delta\phi_{0}^{2}=0.001$. During the first
$\sim 200$ ms, we observe a monotonic increase in the spatial variance, which
can be described by a power law
$\langle\Delta\phi(t)^{2}\rangle-\delta\phi_{0}^{2}\simeq t^{2}$ as shown in
Fig.2(e). This initial growth can be attributed to the dynamics of the
traveling composition layers D0SM00964D . Later the nonlinear growth of a
water-rich droplet around the cap dominates, leading to a complex non-
monotonic time dependence of $\langle\Delta\phi(t)^{2}\rangle$, with a maximum
followed by a subsequent decay until it saturates to a mean constant value
corresponding to the final nonequilibrium steady concentration field. It is
worth noting that this behavior of the spatial variance of $\phi$ is similar
to the initial overshoot in response to an external field found in diverse
soft materials, such as actin networks levin2020kinetics , colloidal
suspensions zia2013stress ; sentjabrskaja2014transient , polymer melts
tseng2021constitutive , or a yield-stress fluids benzi2021stress . In our
system, the maximum in $\langle\Delta\phi(t)^{2}\rangle$ results from local
phase separation around the warmer part of the colloid after switching on
illumination, which proceeds via a two-step process under a step-like
temperature variation across the critical temperature- This mechanism shares
some similarities with the origin of the overshoot in the nonlinear response
of entangled polymer solutions ravindranath2008universal , where the stress
overshoot is also a consequence of a two-step process under a step-like shear
rate. In that case, there is an initial elastic energy storage mechanism
(similar to the transient layering in our system) followed by a sudden
dissipative energy release after the stress peak (in our system, the droplet
formation). Significant fluctuations around the mean steady state value of
$\langle\Delta\phi(t)^{2}\rangle$ occur due to the strong concentration
fluctuations taking place inside the water-rich droplet, as will be discussed
later. The duration of the transient state until the concentration field
relaxes to the nonequilibrium steady state strongly depends on the value of
the $\Delta T$. Interestingly, a further increase in $\Delta T$ gives rise to
a qualitatively different behavior of the time evolution of
$\langle\Delta\phi(t)^{2}\rangle$, in which the transient exhibits a very slow
relaxation to the nonequilibrium steady state with strong concentration
fluctuations and more than a single maximum - see Fig.2(c) and (d). This is
correlated with the final shape of the water-rich droplet formed in the
nonequilibrium steady state. In the case shown in Fig.2(a) and (b), the
droplet covers partially the colloid surface (see the insets). However, for
$\Delta T>9.1^{\circ}$C, the temperature is sufficiently high to form a
droplet that completely surrounds the colloidal particle.

Figure 3: (a) Dependence on the rescaled time
$(t/\tau_{\mathrm{off}})^{\beta}$ with $\beta=0.62$ of the difference between
the spatial variance of $\phi(r,t)$ and the offset value $\delta\phi_{0}^{2}$,
normalized by the amplitude $A$, after switching off the green illumination at
time $t=0$ s, for different values of $\Delta T$ (symbols). The symbols
correspond to: $\Delta T$ = 7.8∘ C (blue circles), 9.1∘ C (green squares),
10.4∘ C (orange diamonds) and 11.7∘ C (red triangles). The solid line
represents Eq.(3). (b) Relaxation time $\tau_{\mathrm{off}}$ of nonequilibrium
concentration profile after switching off the laser for different values of
the temperature difference $\Delta T$.

The relaxation from the nonequilibrium steady concentration field under green
illumination to a thermal-equilibrium state at temperature $T_{0}=28^{\circ}$C
is achieved by switching off the laser at time $t=0$, as shown in Fig. 2. In
this case, the time evolution of $\langle\Delta\phi(t)^{2}\rangle$ is
monotonically decreasing in time for all values of $\Delta T$. In the long
time limit, the spatial variance of the order parameter $\phi(r,t)$ reaches
the constant value $\langle\Delta\phi(t)^{2}\rangle=\delta\phi_{0}^{2}=0.001$,
thus revealing the full relaxation of the fluid to thermal equilibrium, where
the concentration profile is uniform. As shown in Fig.3(a), the relaxation is
well described by the stretched exponential

$\langle\Delta\phi(t)^{2}\rangle=\delta\phi_{0}^{2}+A\exp\left[-\left(\frac{t}{\tau_{\mathrm{off}}}\right)^{\beta}\right],$
(3)

with $\beta=0.62$ for all values of $\Delta T$, regardless of the shape of the
nonequilibrium stationary droplet. As shown in Fig.3(b), the resulting values
of the relaxation time $\tau_{\mathrm{off}}$ decrease upon increasing the
temperature difference $\Delta T$. The stretch exponential, called in physics
Kohlrausch–Williams–Watts (KWW) function wu2016heterogeneous ;
lukichev2019physical , is known to describe relaxations if, e.g., each
constituent particles relaxing exponentially through a single energy barrier,
but there are different barrier heights for different particles such as in
glassy systems xia2001microscopic ; welch2013dynamics , buckled colloidal
monolayers han2008geometric , cytoskeletal networks lieleg2011slow , or the
luminiscence decay of colloidal quantum dots bodunov2017origin . In our
system, such a distribution of different barrier heights for the relaxation of
the concentration field is produced by the spatial heterogeneity created by
initial temperature gradient.

##  III Theoretical description

Our theoretical description of the relaxation process is based on the fluid
particle dynamics method PhysRevLett.85.1338 . In this approach, the Janus
particle is described by the shape function
$\mathcal{S}(\mathbf{r},\mathbf{x})=\frac{1}{2}\left[1+\tanh\left(\frac{R-|\mathbf{r}-\mathbf{x}|}{d_{\mathcal{S}}}\right)\right],$
and the particle orientation function
$\mathcal{O}(\mathbf{r},\mathbf{x},\mathbf{n})=\frac{1}{2}\left[1+\tanh\left(\frac{1}{d_{\mathcal{O}}}\mathbf{n}\cdot\frac{\mathbf{r}-\mathbf{x}}{|\mathbf{r}-\mathbf{x}|}\right)\right],$
where $\mathbf{r}$ is the coordinate of the lattice space and $\mathbf{x}$ is
the position of the center of the particle. $d_{\mathcal{S}}$ represents the
width of the smooth interface such that, in the limit of
$d_{\mathcal{S}}\rightarrow 0$, $\mathcal{S}$ is unity and zero in the
interior and exterior of the particle, respectively. $\mathbf{n}$ is the unit
vector along the particle orientation, and $d_{\mathcal{O}}$ is a sharpness
parameter of the particle orientation. Roughly, one has $\mathcal{O}=1$ on the
capped side
($\mathbf{n}\cdot(\mathbf{r}-\mathbf{x})/|\mathbf{r}-\mathbf{x}|>d_{\mathcal{O}}$).
Otherwise, one has $\mathcal{O}=0$. We consider a binary liquid mixture with
an upper critical temperature of demixing (UCP), therefore heating by
illumination in the experiment corresponds to cooling in our model. In our
modeling, we focus on universal features of the system near the UCP, which are
the same as for the LCP, in order to be able to compare with the experiment.
If one would like to model the LCP beyond the universal features, one would
have to take into account specific interactions, e.g., hydrogen bonding, that
give rise to mixing in low temperatures. The mixture is characterized by the
concentration field $\varphi(\mathbf{r})$, and the temperature field
$T(\mathbf{r})$. The local energy density of the binary mixture
$e(\varphi,\mathbf{x},\mathbf{n},T)$ is a sum of the kinetic energy
$e_{T}(T)=\frac{3}{2}k_{\mathrm{B}}T$ and the interaction energy
$e_{\varphi}(\varphi,\mathbf{x},\mathbf{n})=(1-\mathcal{S})\left[-\frac{\epsilon}{2}\varphi^{2}+\frac{C}{2}(\nabla\varphi)^{2}\right]+\frac{\chi_{\mathrm{p}}}{2}(\varphi-\varphi_{\mathrm{p}})^{2}\mathcal{S}+\left[W_{\mathrm{u}}+(W_{\mathrm{c}}-W_{\mathrm{u}})\mathcal{O}\right]\varphi|d_{\mathcal{S}}\nabla\mathcal{S}|$.
$k_{\mathrm{B}}(=1)$ is the Boltzmann constant and $\epsilon(>0)$ is the
interaction parameter. $W_{\mathrm{c}}$ and $W_{\mathrm{u}}$ represent the
symmetry breaking surface fields. The second term in $e_{\varphi}$ is
introduced to avoid solvent invasion into the particle. $\chi_{\mathrm{p}}$
and $\varphi_{\mathrm{p}}$ are its control parameters. Since $\varphi$ is
conserved locally, its time development is given by the conservation law:

$\displaystyle\frac{\partial\varphi}{\partial t}$ $\displaystyle=$
$\displaystyle-\nabla\cdot(\varphi\mathbf{v})-\nabla\cdot\left(-L_{\varphi}\frac{\mu}{T}\right)+\zeta(\mathbf{r},t),$
(4) $\displaystyle\frac{\partial e}{\partial t}$ $\displaystyle=$
$\displaystyle-\nabla\cdot(e\mathbf{v})-\nabla\cdot\left(L_{T}\nabla\frac{1}{T}\right)-h\mathcal{O}|d_{\mathcal{S}}\nabla\mathcal{S}|,$
(5)

where $L_{\varphi}$ and $L_{T}$ are positive kinetic coefficients. The
Gaussian white noise $\zeta$ obeys the relation
$\langle\zeta(\vec{r},t)\leavevmode\nobreak\
\zeta(\vec{r}^{\prime},t^{\prime})\rangle=-2\zeta_{0}\nabla^{2}\delta(\vec{r}-\vec{r}^{\prime})\delta(t-t^{\prime})$;
$\zeta_{0}$ is the strength of the noise. We ignore the off-diagonal terms.
$h$ represents the cooling power, and the temperature field is calculated from
$e(\varphi,\mathbf{x},\mathbf{n},T)$. Because the system is nonisothermal, the
chemical potential is calculated from the entropy density $s$ as
$\mu=-T\left(\frac{\partial s}{\partial\varphi}\right)_{e}$. Expressions for
$s$ and $\mu$ are given in Appendix A. In mean field, $\mu$ gives the critical
point at $\varphi_{\mathrm{c}}=0$ and
$T_{\mathrm{c}}=\epsilon/(ak_{\mathrm{B}})$. When $T<T_{\mathrm{c}}$, the UCP
symmetric mixture is phase-separated. The velocity $\mathbf{v}(\mathbf{r})$ of
the hydrodynamic flow obeys the following hydrodynamic equation:

$\displaystyle
C\nabla\cdot(\nabla\varphi:\nabla\varphi)-\frac{\mathcal{S}}{\Omega}\frac{\partial
e}{\partial\mathrm{x}}-\frac{1}{2}\nabla\times\left\\{\frac{\mathcal{S}}{\Omega}\mathbf{n}\times\left(\frac{\partial
e}{\partial\mathbf{n}}\right)\right\\}-\nabla p$
$\displaystyle+\nabla\left[\\{\eta+(\eta_{\mathrm{c}}-\eta)\mathcal{S}\\}\\{\nabla:\mathbf{v}+(\nabla:\mathbf{v})^{T}\right]+\mathbf{F_{x}}=0.$
(6)

The first term is the mechanical stress stemming from the concentration
inhomogeneity, i.e., the interface tension. The second and third terms are due
to the particle translation and rotation. $p$ is the pressure obtained by the
incompressible condition $\nabla\cdot\mathbf{v}=0$. Within the FPD method, the
fifth term is due to the viscous stress, in which $\eta$ and
$\eta_{\mathrm{c}}$ are the viscosity of the solvent and the inside of the
particle, respectively. The last term
$\mathbf{F_{x}}=K(\mathbf{x}-\mathbf{x}_{0})\frac{\mathcal{S}}{\Omega}$ is
introduced in order to fix the particle at its initial position
$\mathbf{x}_{0}$ by imposing a harmonic potential with spring constant $K$. By
neglecting the convection terms in Eqs.(4) and (5) the model reduces to the
nonisothermal purely diffusive model of type B.

Figure 4: Numerical results of the hydrodynamic model: (a) Snapshot of the
concentration field around the Janus particle at time $t=-400$. The capped
surface (oriented toward the $x$ axis) is cooled with the power $h=4$; (b)
Time evolution of the spatial variance $\langle\Delta\varphi(t)^{2}\rangle$ of
the concentration field averaged over the whole sample for $h=4$; (c) Time
evolution of the fluctuations $\langle v^{2}\rangle$ of the velocity field
averaged over the whole sample for several values of the cooling power $h$.

##  IV Numerical results

We solve the evolution equations numerically (for details of numerical
simulation see Appendix B). For the initial configuration we assume
$\varphi(r,t)=0$ throughout, $T=1.1T_{\mathrm{c}}$ everywhere and no flow. At
$t=-400$ we perform a temperature quench at the capped surface of the particle
by switching on the cooling power $h$. At $t=0$ the cooling power is switched
off and the time evolution is registered up to $t=400$. The snapshot of the
concentration field $\varphi(r,t)$ at $t=0$ and $h=4$ is shown in Fig.4(a).
The binary solvent is phase separated around the colloid forming an asymmetric
droplet within which the concentration field strongly fluctuates. Fig.4(b)
shows how the spatial variance $\langle\Delta\varphi(t)^{2}\rangle$ of the
concentration field averaged over the whole sample changes in the whole period
of time for various temperature quenches. Velocity fluctuations are shown in
Fig.$\ref{Fig:Hprofiles}(c)$.

Figure 5: Numerical results of the hydrodynamic model: (a) Fitting numerical
results for normalized spatial variance
$\langle\Delta\varphi(r,t)(t)^{2}\rangle_{1D}$ of the concentration field
calculated along the main particle axis during the OFF process to the
stretched exponential decay given by Eq.(3) (dashed line) with $\beta=0.62$,
for different values of $h$ (symbols). $\delta\varphi_{0}=0.00303$. (b)
Relaxation time $\tau_{\mathrm{off}}$ of nonequilibrium concentration profile
after switching off cooling for different values of the cooling power $h$.

During the ON process $\langle\Delta\varphi(t)^{2}\rangle$ resulting from the
hydrodynamic model exhibits a noisy nonmonotonic behavior which is very
similar to the behavior of $\langle\Delta\phi(t)^{2}\rangle$ observed in
experiment. On the contrary, purely diffusive model B gives a smooth and
monotonic evolution towards the steady state (see Fig.7(a) in Appendix C). It
is well established that thermal fluctuations in liquids in the presence of
stationary temperature gradients are anomalously large and very long-ranged
de2006hydrodynamic . They occur as a result of a coupling between temperature
and velocity fluctuations. In liquid mixtures a stationary temperature
gradient induces stationary concentration gradient via the Soret effect.
Fluctuating hydrodynamics for a binary fluid mixtures bounded between two
parallel plates with different temperatures predicts that nonequilibrium
fluctuations of concentration, temperature, as well as velocity will be
present at a stationary state of a mixture croccolo2016non . However, because
in liquid mixtures the ratio between thermal diffusivity and the mass
diffusion coefficient (Lewis number) and the ratio between the kinematic
viscosity and the mass diffusion coefficient (Schmidt number) are commonly
larger than unity, nonequilibrium fluctuations in concentration will be
dominant. Here we study nonequilibrium fluctuations not only at a stationary
state, but also when approaching it. In the initial stage, the hydrodynamic
flow is weak so that both hydrodynamic and purely diffusive models give a
similar growth of concentration fluctuations:
$\langle\Delta\varphi(r,t)(t)^{2}\rangle_{1D}$ calculated along the main
particle axis $R<r<R+L$ with $L=100$ in the $x$–$y$ plane passing through the
particle can be described by the power law $\sim t^{2}$ (see Figs.2(f) and
7(b) in Appendix C), in a perfect agreement with the experimental data. We
stress that this initial behavior cannot be described by linearized equations.
In the late stage, on the other hand, the concentration fluctuations are
caused by coupling between the concentration and velocity fields, because the
Lewis number is large in our simulations. From Figs.4(b) and (c), we can see
that for our choice of parameters the velocity fluctuations are comparable
with the concentration fluctuations. Other predictions of fluctuating
hydrodynamics indicate that the intensity of nonequilibrium fluctuations is
proportional to $\left(\nabla T\right)$ and that for small wave numbers $q$
the intensity of the nonequilibrium concentration fluctuations diverges as
$q^{-4}$, which is much stronger than the divergence of critical concentration
fluctuations as $q^{-2}$ near a consolute point croccolo2016non . Comparing
our experimental data shown in Figs.2(a) and (d) and the theoretical curves
for $h=1$ and $h=4$ in Figs.4(b) and (c) we can conjecture that indeed,
thermal fluctuations are stronger for larger temperature gradients. Due to
technical difficulties, we have not analyzed the $q$-dependence of our
results.

Figure 6:  Evolution of (a) the concentration field contribution $\langle
e_{\varphi}\rangle$ to the energy and of (b) the total energy $\langle
e\rangle$ of the system for different values of the cooling power $h$.
$\langle e\rangle$ and $\langle e_{\varphi}\rangle$ are the spatial averages
of $e$ and $e_{\varphi}$.(c) Initial decay time $\tau_{0}$ of the total energy
$\langle e\rangle$ after switching off the laser and the ratio $\Delta\langle
e(t=0)\rangle V/J$ between the energy stored in the steady state and the total
energy flux at this state for different values of $h$. $V$ is the system
volume.

In the simulations, we also observe that returning to equilibrium is much
faster than approaching the steady state. Remarkably, during the OFF process
the time dependence of $\langle\varphi(r,t)(t)^{2}\rangle_{1D}$ is well
described by Eq.(3) with the same exponent $\beta=0.62$ as in the experiment,
as shown in Fig.5. The characteristic time of the stretched exponential decay
decreases with the cooling power. This trend is also similar to the one
observed in the experiment. With increasing values of $h$ or $\Delta T$, a
larger temperature gradient builds up around the colloid, which speeds up
relaxation to the equilibrium. Next, we compute the energy of the system,
which is practically not accessible from experiment. As can be inferred from
Fig.6(a) and (b), $\langle e\rangle$ is dominated by the kinetic contribution
$\langle e_{T}\rangle$ and hence it decreases after the cooling is switched on
and increases upon switching off the cooling power. Restructuring of the
concentration field in response to the temperature changes is much slower in
both ON and OFF processes, as mirrored in the evolution of $\langle
e_{\varphi}\rangle$. Achieving the nonequilibrium stationary value of $\langle
e_{\varphi}\rangle$ takes longer because the initial concentration field is
uniformly distributed whereas returning to equilibrium starts from a
concentrated distribution and therefore is faster (see movies in SM). Finally,
relaxation of the energy to equilibrium is well described by Eq.(3) with the
exponent $\beta=0.62$ and the characteristic decay time $\tau_{\textrm{off}}$
which decreases with $h$ as can be seen in Fig.6(c). Decay of
$\langle\varphi(r,t)(t)^{2}\rangle_{1D}$ and $\langle e_{T}\rangle$ in the
purely diffusive model B is better described by a sum of the exponential end
stretched exponential functions of time - see Appendix C. We find that
$\tau_{\textrm{off}}$ is almost the same as the initial decay time $\tau_{0}$
obtained from the change rate of the energy at $t=0$. Interestingly, both
characteristic times are proportional to the ratio $\mathcal{T}$ between the
energy stored in the nonequilibrium steady state and the total energy flux in
this state $\mathcal{T}=\Delta\langle e(t=0)\rangle V/J$\- see Fig.6(c). The
same holds for the initial decay time $\tau_{0}$ obtained in model B, as can
be seen in Fig.8(d) in Appendix C. Similar observation was made for the energy
out-flow from the Lennard-Jones system - but only for initial times after the
shutdown of energy flux into the system PhysRevE.99.042118 .

##  V Summary

In summary, we have investigated the relaxation process of the concentration
field of a binary mixture around a Janus colloid heated by laser illumination
through the critical temperature of the fluid. Our results show a pronounced
asymmetry in the approach to a nonequilibrium state under a temperature
gradient around the colloid with respect to the relaxation of the fluid back
to equilibrium at constant temperature. While in the former the dynamics is
dominated by transient composition layers until the final noisy steady state
is reached, the latter is faster and follows a stretched-exponential decay
consistent with the characteristic timescale of energy outflow after the
shutdown of the energy input. In particular, we show that such relaxation
timescale is directly related to the energy stored in the nonequilibrium
steady state and the total energy flux. Further experimental and theoretical
efforts could help to address the connection between these two quantities with
the relaxation from nonequilibrium states in other systems.

##  VI Acknowledgments

A. M. thanks Sutapa Roy for discussions on the purely diffusive model. T.A. is
financially supported by KAKENHI (Grant No.21K03486), and CREST, JST
(JPMJCR2095). J. R. G.-S. was supported by UNAM-PAPIIT IA103320

##  Appendix A Entropy density and chemical potential of the hydrodynamic
model

Close to the critical point the entropy density can be written as

$s(\varphi,\mathbf{x},T)=-k_{\mathrm{B}}(1-\mathcal{S})\left(\frac{a}{2}\varphi^{2}+\frac{b}{4}\varphi^{4}\right)+k_{\mathrm{B}}\left[\ln\left(\frac{v}{\lambda_{T}^{3}}\right)+1\right],$
(7)

where $a$ and $b$ are positive constants, $\lambda_{T}(\propto T^{-1/2}$) is
the thermal de Broglie length, and $v$ is the molecular volume, thus

$\displaystyle\mu$ $\displaystyle=$
$\displaystyle(1-\mathcal{S})\left[(k_{\mathrm{B}}Ta-\epsilon)\varphi+k_{\mathrm{B}}Tb\varphi^{3}\right]-C\nabla\cdot\left[(1-\mathcal{S})\nabla\varphi\right]$
(8)
$\displaystyle+[W_{\mathrm{u}}+(W_{\mathrm{c}}-W_{\mathrm{u}})\mathcal{O}]|d_{\mathcal{S}}\nabla\mathcal{S}|+\chi_{\mathrm{p}}(\varphi-\varphi_{\mathrm{p}})\mathcal{S}.$

##  Appendix B Parameters of the numerical simulations

The simulation box has size $240\times 240\times 52$ in units of the lattice
spacing $\Delta=1$. The walls are placed at $z=1$ and $z=52$. The non-slip
boundary conditions for the hydrodynamic flow are imposed at the walls, i.e.,
$\mathbf{v}(z=1)=\mathbf{v}(z=52)=0$. The temperatures at the walls are kept
above the critical temperature $T(z=1)=T(z=52)=1.1T_{\mathrm{c}}$. The
concentration flux $-L_{\varphi}\nabla(\mu/T)$ is imposed to vanish at the
walls. The periodic boundary conditions are applied to the $x$ and $y$
directions. The time increment is $\Delta t=0.001$. We take the particle
radius $R=20$ and set $d_{\mathcal{S}}=1$ and $d_{\mathcal{O}}=0.033$. The
capped surface is strongly hydrophilic with $W_{\mathrm{c}}=-2.0$, while the
uncapped surface is weakly hydrophilic with $W_{\mathrm{u}}=-0.2$. In Eq. (8)
we set $\epsilon=1$, $C=4$, $\chi_{\mathrm{p}}=5$, $\varphi_{\mathrm{p}}=0$,
$a=1$ and $b=1$. In equations for time evolution (Eqs. (4) and (5)), we employ
$L_{\varphi}=(1-\mathcal{S})$, $L_{T}=10(1+2\mathcal{S})$, and
$\zeta_{0}=10^{-4}$. In the hydrodynamic equation (6), we set $\eta=0.5$,
$\eta_{\mathrm{c}}=25$ and $K=100$.

Calculation of stationary state fluxes The inlet and outlet fluxes are
calculated as

$\displaystyle J_{\mathrm{in}}$ $\displaystyle=$
$\displaystyle\int\mathrm{d}\mathbf{r}h\mathcal{O}|d_{\mathcal{S}}\nabla\mathcal{S}|,$
(9) $\displaystyle J_{\mathrm{out}}$ $\displaystyle=$
$\displaystyle\int\mathrm{d}\mathbf{A}\cdot L_{T}\nabla\frac{1}{T},$ (10)

where $\int\mathrm{d}\mathbf{A}$ represents the integral over the two walls.

##  Appendix C Results from the purely diffusive model B

Figure 7: Numerical results from the purely diffusive model B. (a) Time
evolution for the spatial variance $\langle\Delta\varphi(t)^{2}\rangle$ of the
concentration field averaged over the whole system for different values of the
cooling power $h$; (b) The initial growth of
$\langle\Delta\varphi(t)^{2}\rangle_{1D}$ of the concentration field
calculated along the main particle axis after switching on cooling at
$t=-5000$. $\delta\varphi_{0}^{2}=0.00303$; (c) Normalized
$\langle\Delta\varphi(t)^{2}\rangle_{1D}$ after switching off cooling at time
$t=0$ for different values of $h$. The horizontal axis is scaled with the
characteristic decay time of the fast mode (see Eq. (11)). (d) The $h$
dependence of two times $\tau_{\textrm{fast}}$ and $\tau_{\textrm{slow}}$
characterizing the decay of $\langle\Delta\varphi(t)^{2}\rangle_{1D}$ after
switching off cooling obtained from fitting to the function $f(t)$ given by
Eq. (11) with $\beta\approx 0.62$.  Figure 8: Numerical results for time
evolution of (a) the contribution $\langle e_{\varphi}\rangle$ to the total
energy density due to the concentration field and of (b) the total energy
density $\langle e\rangle$ of the system calculated within the purely
diffusive model B for different values of the cooling power $h$. $\langle
e\rangle$ and $\langle e_{\varphi}\rangle$ are the spatial averages of $e$ and
$e_{\varphi}$, where $V$ is the system volume. (c) The $h$ dependence of two
times $\tau_{\textrm{fast}}$ and $\tau_{\textrm{slow}}$ characterizing the
decay of $\langle e\rangle$ after switching off cooling obtained from fitting
to the function $f(t)$ given by Eq. (11). (d) Initial decay time $\tau_{0}$ of
$\langle e\rangle$ after switching off the laser and the ratio $\Delta\langle
e(t=0)\rangle V/J$ between the energy stored in the steady state and the total
energy flux at this state for different values of $h$.

The cooling power is applied at $t=-5000$ and removed at $t=0$. The spatial
variance $\langle\Delta\varphi(t)^{2}\rangle$ of the concentration field
averaged over the whole system increases monotonically until the cooling is
switched off (Fig. 7(a)). After a rapid initial growth, which is due to
concentration waves traveling from the surface,
$\langle\Delta\varphi(t)^{2}\rangle$ approaches the steady state value
logarithmically. The initial growth of the concentration field
$\langle\Delta\varphi(t)^{2}\rangle_{1D}$ calculated along the main particle
axis follows well the power law $\approx t^{2}$, in agreement with
experimental and hydrodynamic model results (Fig. 7(b)).

The decay of $\langle\Delta\varphi(t)^{2}\rangle_{1D}$ after the removal of
$h$ is fitted well by the function

$f(t)=a_{1}\exp(-t/\tau_{\textrm{fast}})+a_{2}\exp(-(t/\tau_{\textrm{slow}})^{\beta})+c$
(11)

with experimental value $\beta\approx 0.62$. Characteristic time
$\tau_{\textrm{fast}}$ of fast exponential decay dominating at very short
times is almost independent of the cooling power. Time $\tau_{\textrm{slow}}$
characterizing the stretched exponential decay decreases with increasing
cooling power. If we treat $\beta$ in Eq. (11) as a fitting parameter, both
characteristic times change only slightly. The fitted values of $\beta$ varies
from $\approx 0.72$ for $h=0.5$ to $\approx 0.62$ for $h=4$.

Results for the total energy density $\langle e\rangle$ averaged over the
system and the contribution to it from the concentration field $\langle
e_{\varphi}\rangle$ are shown in Fig. 8. After switching off cooling $e$ is
also well described by Eq. (11) with the stretched exponent $\beta\approx$
0.58.

## References

  * (1) L. Onsager, “Reciprocal relations in irreversible processes. i.”, Phys. Rev. 37, 405 (1931).
  * (2) L. Onsager, “Reciprocal relations in irreversible processes. ii.”, Phys. Rev. 38, 2265 (1931).
  * (3) R. Kubo, M. Yokota, and S. Nakajima, “Statistical-mechanical theory of irreversible processes. ii. response to thermal disturbance”, J. Phys. Soc 12, 1203 (1957).
  * (4) S. Dattagupta, Relaxation phenomena in condensed matter physics, Elsevier (2012).
  * (5) R. Metzler, E. Barkai, and J. Klafter, “Anomalous diffusion and relaxation close to thermal equilibrium: A fractional fokker-planck equation approach”, Phys. Rev. Lett. 82, 3563 (1999).
  * (6) N. Shiraishi and K. Saito, “Information-theoretical bound of the irreversibility in thermal relaxation processes”, Phys. Rev. Lett. 123, 110603 (2019).
  * (7) A. Lapolla and A. Godec, “Faster uphill relaxation in thermodynamically equidistant temperature quenches”, Phys. Rev. Lett. 125, 110602 (2020).
  * (8) C. Maes, K. Netočný, and B. Wynants, “Monotonic return to steady nonequilibrium”, Phys. Rev. Lett. 107, 010601 (2011).
  * (9) L. V. F. Gomes and A. B. Kolomeisky, “Dynamics of relaxation to a stationary state for interacting molecular motors”, J. Phys. A 51, 015601 (2017).
  * (10) A. Dhar, A. Kundu, S. N. Majumdar, S. Sabhapandit, and G. Schehr, “Run-and-tumble particle in one-dimensional confining potentials: Steady-state, relaxation, and first-passage properties”, Phys. Rev. E 99, 032132 (2019).
  * (11) J. R. Gomez-Solano, A. Petrosyan, S. Ciliberto, R. Chetrite, and K. Gawędzki, “Experimental verification of a modified fluctuation-dissipation relation for a micron-sized particle in a nonequilibrium steady state”, Phys. Rev. Lett. 103, 040601 (2009).
  * (12) V. Blickle, J. Mehl, and C. Bechinger, “Relaxation of a colloidal particle into a nonequilibrium steady state”, Phys. Rev. E 79, 060104 (2009).
  * (13) J. Gieseler, L. Novotny, C. Moritz, and C. Dellago, “Non-equilibrium steady state of a driven levitated particle with feedback cooling”, New J. Phys. 17, 045011 (2015).
  * (14) J. A. Owen, T. R. Gingrich, and J. M. Horowitz, “Universal thermodynamic bounds on nonequilibrium response with biochemical applications”, Phys. Rev. X 10, 011066 (2020).
  * (15) U. Basu, L. Helden, and M. Krüger, “Extrapolation to nonequilibrium from coarse-grained response theory”, Phys. Rev. Lett. 120, 180604 (2018).
  * (16) G. Baffou and R. Quidant, “Thermo-plasmonics: using metallic nanostructures as nano-sources of heat”, Laser Photon. Rev. 7, 171 (2013).
  * (17) A. Torchi, F. Simonelli, R. Ferrando, and G. Rossi, “Enhancement of lipid membrane permeability induced by irradiated gold nanoparticles”, ACS Nano 11, 12553 (2017).
  * (18) A. Bahadori, G. Moreno-Pescador, L. B. Oddershede, and P. M. Bendix, “Remotely controlled fusion of selected vesicles and living cells: a key issue review”, Rep. Prog. Phys. 81, 032602 (2018).
  * (19) G. Volpe, I. Buttinoni, D. Vogt, H.-J. Kümmerer, and C. Bechinger, “Microswimmers in patterned environments”, Soft Matter 7, 8810 (2011).
  * (20) I. Buttinoni, G. Volpe, F. Kümmel, G. Volpe, and C. Bechinger, “Active brownian motion tunable by light”, J. Phys.: Condens. Matter 24, 284129 (2012).
  * (21) J. R. Gomez-Solano, A. Blokhuis, and C. Bechinger, “Dynamics of self-propelled janus particles in viscoelastic fluids”, Phys. Rev. Lett. 116, 138301 (2016).
  * (22) S. Kumar, A. Kumar, M. Gunaseelan, R. Vaippully, D. Chakraborty, J. Senthilselvan, and B. Roy, “Trapped in out-of-equilibrium stationary state: Hot brownian motion in optically trapped upconverting nanoparticles”, Front. Phys. 8, 429 (2020).
  * (23) H. R. Vutukuri, M. Lisicki, E. Lauga, and J. Vermant, “Light-switchable propulsion of active particles with reversible interactions”, Nat. Commun. 11, 1 (2020).
  * (24) S. Auschra, A. Bregulla, K. Kroy, and F. Cichos, “Thermotaxis of janus particles”, Eur. Phys. J. E 44, 90 (2021).
  * (25) A. Girot, N. Danné, A. Würger, T. Bickel, F. Ren, J. Loudet, and B. Pouligny, “Motion of optically heated spheres at the water–air interface”, Langmuir 32, 2687 (2016).
  * (26) F. Schmidt, A. Magazzù, A. Callegari, L. Biancofiore, F. Cichos, and G. Volpe, “Microscopic engine powered by critical demixing”, Phys. Rev. Lett. 120, 068004 (2018).
  * (27) F. Schmidt, H. Šípová-Jungová, M. Käll, A. Würger, and G. Volpe, “Non-equilibrium properties of an active nanoparticle in a harmonic potential”, Nat. Commun. 12, 1 (2021).
  * (28) A. Onuki, Phase transition dynamics, Cambridge University Press (2002).
  * (29) M. Shono, S. Takatori, J. M. Carnerero, and K. Yoshikawa, “Switching between positive and negative movement near an air/water interface through lateral laser illumination”, Appl. Phys. Lett. 117, 073701 (2020).
  * (30) M. Pascual, A. Poquet, A. Vilquin, and M.-C. Jullien, “Phase separation of an ionic liquid mixture assisted by a temperature gradient”, Phys. Rev. Fluids 6, 024001 (2021).
  * (31) T. Zinn, L. Sharpnack, and T. Narayanan, “Phoretic dynamics of colloids in a phase separating critical liquid mixture”, Phys. Rev. Research 2, 033177 (2020).
  * (32) C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, “Active particles in complex and crowded environments”, Rev. Mod. Phys. 88, 045006 (2016).
  * (33) S. Roy, S. Dietrich, and A. Maciołek, “Solvent coarsening around colloids driven by temperature gradients”, Phys. Rev. E 97, 042603 (2018).
  * (34) S. Roy and A. Maciołek, “Phase separation around a heated colloid in bulk and under confinement”, Soft Matter 14, 9326 (2018).
  * (35) A. Würger, “Self-diffusiophoresis of janus particles in near-critical mixtures”, Phys. Rev. Lett. 115, 188304 (2015).
  * (36) J. R. Gomez-Solano, S. Samin, C. Lozano, P. Ruedas-Batuecas, R. van Roij, and C. Bechinger, “Tuning the motility and directionality of self-propelled colloids”, Sci. Rep. 7, 1 (2017).
  * (37) T. Araki and A. Maciołek, “Illumination-induced motion of a janus nanoparticle in binary solvents”, Soft Matter 15, 5243 (2019).
  * (38) J. R. Gomez-Solano, S. Roy, T. Araki, S. Dietrich, and A. Maciołek, “Transient coarsening and the motility of optically heated janus colloids in a binary liquid mixture”, Soft Matter 16, 8359 (2020).
  * (39) P. Bauduin, L. Wattebled, S. Schrödle, D. Touraud, and W. Kunz, “Temperature dependence of industrial propylene glycol alkyl ether/water mixtures”, J. Mol. Liq. 115, 23 (2004).
  * (40) M. Assenheimer, B. Khaykovich, and V. Steinberg, “Phase separation of a critical binary mixture subjected to a temperature gradient”, Phys. A: Stat. Mech. Appl. 208, 373 (1994).
  * (41) C. Mauger, L. Méès, M. Michard, A. Azouzi, and S. Valette, “Shadowgraph, schlieren and interferometry in a 2d cavitating channel flow”, Exp. in Fluids 53, 1895 (2012).
  * (42) M. Levin, R. Sorkin, D. Pine, R. Granek, A. Bernheim-Groswasser, and Y. Roichman, “Kinetics of actin networks formation measured by time resolved particle-tracking microrheology”, Soft Matter 16, 7869 (2020).
  * (43) R. N. Zia and J. F. Brady, “Stress development, relaxation, and memory in colloidal dispersions: Transient nonlinear microrheology”, J. Rheol. 57, 457 (2013).
  * (44) T. Sentjabrskaja, M. Hermes, W. Poon, C. Estrada, R. Castaneda-Priego, S. Egelhaaf, and M. Laurati, “Transient dynamics during stress overshoots in binary colloidal glasses”, Soft Matter 10, 6546 (2014).
  * (45) H.-C. Tseng, “A constitutive analysis of stress overshoot for polymer melts under startup shear flow”, Phys. Fluids 33, 051706 (2021).
  * (46) R. Benzi, T. Divoux, C. Barentin, S. Manneville, M. Sbragaglia, and F. Toschi, “Stress overshoots in simple yield stress fluids”, Phys. Rev. Lett. 127, 148003 (2021).
  * (47) S. Ravindranath and S.-Q. Wang, “Universal scaling characteristics of stress overshoot in startup shear of entangled polymer solutions”, J. Rheol. 52, 681 (2008).
  * (48) J. H. Wu and Q. Jia, “The heterogeneous energy landscape expression of kww relaxation”, Sci. Rep 6, 1 (2016).
  * (49) A. Lukichev, “Physical meaning of the stretched exponential kohlrausch function”, Phys. Lett. A 383, 2983 (2019).
  * (50) X. Xia and P. G. Wolynes, “Microscopic theory of heterogeneity and nonexponential relaxations in supercooled liquids”, Phys. Rev. Lett. 86, 5526 (2001).
  * (51) R. C. Welch, J. R. Smith, M. Potuzak, X. Guo, B. F. Bowden, T. Kiczenski, D. C. Allan, E. A. King, A. J. Ellison, and J. C. Mauro, “Dynamics of glass relaxation at room temperature”, Phys. Rev. Lett. 110, 265901 (2013).
  * (52) Y. Han, Y. Shokef, A. M. Alsayed, P. Yunker, T. C. Lubensky, and A. G. Yodh, “Geometric frustration in buckled colloidal monolayers”, Nature 456, 898 (2008).
  * (53) O. Lieleg, J. Kayser, G. Brambilla, L. Cipelletti, and A. R. Bausch, “Slow dynamics and internal stress relaxation in bundled cytoskeletal networks”, Nat. Mater. 10, 236 (2011).
  * (54) E. Bodunov, Y. A. Antonov, and A. L. Simões Gamboa, “On the origin of stretched exponential (kohlrausch) relaxation kinetics in the room temperature luminescence decay of colloidal quantum dots”, J. Chem. Phys. 146, 114102 (2017).
  * (55) H. Tanaka and T. Araki, “Simulation method of colloidal suspensions with hydrodynamic interactions: Fluid particle dynamics”, Phys. Rev. Lett. 85, 1338 (2000).
  * (56) J. M. O. De Zárate and J. V. Sengers, Hydrodynamic fluctuations in fluids and fluid mixtures, Elsevier (2006).
  * (57) F. Croccolo, J. O. de Zárate, and J. Sengers, “Non-local fluctuation phenomena in liquids”, Eur. Phys. J. E 39, 1 (2016).
  * (58) R. Holyst, A. Maciołek, Y. Zhang, M. Litniewski, P. Knychała, M. Kasprzak, and M. Banaszak, “Flux and storage of energy in nonequilibrium stationary states”, Phys. Rev. E 99, 042118 (2019).


# Effects of Shock Waves on Neutrino Oscillations in Three Supernova Models

Jing Xu1111Email<EMAIL_ADDRESS>Li-Jun Hu 1222Email:
<EMAIL_ADDRESS>Rui-Cheng Li1333Email: Rui-cheng<EMAIL_ADDRESS>Xin-Heng
Guo1444Corresponding author, Email<EMAIL_ADDRESS>and Bing-Lin
Young2,3555Email<EMAIL_ADDRESS>1College of Nuclear Science and
Technology, Beijing Normal University, Beijing 100875, China
2Department of Physics and Astronomy, Iowa State University, Ames, Iowa 5001,
USA
3Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing
100190, China

###### Abstract

It has been realized that the shock wave effects play an important role in
neutrino oscillations during the supernova explosion. In recent years, with
the development of simulations about supernova explosion, we have a better
understanding about the density profiles and the shock waves in supernovae
than before. It has been shown that the appearance of shock waves not only
varies with time, but is also affected by the mass of the supernova. When the
mass of the supernova happens to be in a certain range (e.g. it equals 10.8
times the mass of the sun), there might be a reverse shock wave, another
sudden change of density except the forward shock wave, emerging in the
supernova. In addition, there are some other time-dependent changes of density
profiles in different supernova models. Because of these complex density
profiles, the expression of the crossing probability at the high resonance,
$P_{H}$, which we used previously would be no longer applicable. In order to
get more accurate and reasonable results, we use the data of density profiles
in three different supernova models obtained from simulations to study the
variations of $P_{s}$ (the survival probability of
$\nu_{e}\rightarrow\nu_{e}$), as well as $P_{c}$ (the conversion probability
of $\nu_{x}\rightarrow\nu_{e}$). It is found that the mass of the supernova
does make a difference on the behavior of $P_{s}$, and affects $P_{c}$ at the
same time. With the results of $P_{s}$ and $P_{c}$, we can estimate the number
of $\nu_{e}$ remained after they go through the matter in the supernova.

###### pacs:

14.60.Pq, 13.15.+g, 25.30.Pt, 26.30.-k

##  I Introduction

Since 1980’s, supernova neutrinos have been a focus of our attention for a
long time Arnett1 -Bethe1 . In recent years, the issues about neutrino
detection experiments and simulations of supernova explosion have always been
hotspots in scientific research fields Ardellier -Janka0 . Thanks to the
results of simulations and development of theories on supernovae, the whole
process of supernova explosion has been understood much better now Bruenn
-Kotake . The post bounce period, one stage in the process of supernova
explosion, can be divided generally into the accretion phase and the cooling
phase. During the supernova explosion, both the collective effects and the
shock wave effects play important roles in neutrino oscillations Hudepohl
Chakraborty0 . As we know, in the post bounce period, a large number of
supernova neutrinos go through the supernova matter, carrying an enormous
amount of energy, and emit from the supernova Spergel Giunti . Meanwhile, the
spreading outward energy makes significant changes to the supernova matter and
arouses the shock wave, which has an impact on neutrino oscillations in the
supernova Fogli1 -Schirato .

The essence of the shock wave effect is the Mikheyev-Smirnov-Wolfenstein (MSW)
effect, caused by the interactions of neutrinos with matter and dominated by
the density profile of the supernova as well as neutrino mixing angles
Takahashi1 -Kuo . The neutrino crossing probability in the high resonance,
$P_{H}$, will experience complex changes resulting from the high-density and
unstable matter in the supernova Kuo -Kachelriess . In the past, the
simplified model of supernova density profile, which only contains the forward
shock wave, was used to study the behavior of $P_{H}$ Fogli2 . With this
method, we could roughly figure out the impact of the shock wave effects on
neutrino oscillations in the supernova Guo -Xu . Afterwards, it was found that
besides the forward shock wave there would be a reverse shock wave and some
turbulence existing in the supernova during the process of explosion, all of
which have influences on neutrino oscillations Kneller -Couch . Presently,
from the results of simulations aiming at different supernova models, it is
known that the actual density profiles in supernova are distinctly different
from one to another due to the difference of their masses Lund Wongwathanarat
. For this reason, the original simplified model of the density profile cannot
provide the images of time-varying densities in different supernova models
comprehensively. Therefore, one has to find out the images in other ways. One
possible way is to get new simplified density profiles according to the
simulation results of different supernova models, respectively. In this way,
it is necessary to do data analysis and fitting for densities, which is beyond
the scope of the present work. Therefore, we decide to use the simulation
results of supernova density profiles to discuss the variations of $P_{s}$
(the survival probability of $\nu_{e}\rightarrow\nu_{e}$) and $P_{c}$ (the
conversion probability of $\nu_{x}\rightarrow\nu_{e}$, $\nu_{x}$ represents
$\nu_{\mu}$ and $\nu_{\tau}$) instead of applying the simplified model of
density profiles to study the behavior of $P_{H}$.

This paper is divided into five parts. In Section 2, we give a brief review
about neutrino oscillations in matter. In Section 3, we use the simulation
results of three supernova models to give the images of time-varying density
profiles in order to find the similarities and differences among them. In
Section 4, based on the theory in Section 2 and using the density data in
Section 3, we calculate $P_{s}$, $P_{c}$ and the ratio of $\nu_{e}$ remained
after neutrinos go through the matter in three supernova models. In addition,
we compare the diverse appearances of shock waves among these supernova models
and extract how the time-dependent matter affects neutrino oscillations in
supernovae. The last section consists of a summery and some discussions of
this paper.

##  II Neutrino oscillations in matter

Neutrinos might be regarded as either a Dirac or a Majorana particle. No
matter which kind of particle it is, eliminating the spin structure from its
propagation equation yields the Klein-Gordon equation Kuo . In the flavor
basis, different from the case in vacuum, the propagation equation of neutrino
in matter has a potential energy item, which stems from the interactions
between neutrinos and matter,

$\displaystyle-i\frac{d}{dt}\left|{{\nu}_{\alpha}}\right\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{2p}M^{2}\left|{{\nu}_{\alpha}}\right\rangle,$ (1)
$\displaystyle M^{2}$ $\displaystyle=$
$\displaystyle\left[U\left[\begin{array}[]{cc}m_{1}^{2}&0\\\ 0&m_{2}^{2}\\\
\end{array}\right]U^{\dagger}+\left[\begin{array}[]{cc}A&0\\\ 0&0\\\
\end{array}\right]\right].$ (6)

In Eq. (1), ${\nu}_{\alpha}$ stands for neutrinos in flavor eigenstate, $U$
the unitary transformation between the flavor and mass bases, $m_{1}^{2}$ and
$m_{2}^{2}$ the square of mass eigenvalues, and $A$ is the parameter
containing the effect of matter on neutrino propagation, which is give by

$\displaystyle A=2Vp=2\sqrt{2}G_{F}N_{e}p,$ (7)

where $G_{F}$ is the Fermi constant, $N_{e}$ the number of electrons, and $p$
the momentum of neutrino. For antineutrinos, $A{\rightarrow}-A$ and
$U{\rightarrow}U^{*}$.

For neutrinos, $p\sim E$, so the expression of $A$ can be written as

$\displaystyle A=2\sqrt{2}G_{F}N_{e}E=2\sqrt{2}G_{F}(Y_{e}/m_{n})\rho E,$ (8)

where $\rho$ is the matter density in the supernova, $Y_{e}$ the number of
electron per nucleon, and $m_{n}$ the nucleon mass.

In the study of supernova neutrinos, we work in the approximation that there
are only two neutrino flavors, $\nu_{e}$ and $\nu_{x}$ ($\nu_{x}$ represents
$\nu_{\mu}$ and $\nu_{\tau}$) Kuo . Acting like an induced mass squared of the
electron neutrino, $A$ arises from the propagation of the electron neutrino
through a background of electrons. With $c=1$ and $t\rightarrow x$, the
equation of two flavors in matter is written as

$\displaystyle i\frac{d}{dx}\left[\begin{array}[]{c}{\nu}_{e}\\\ {\nu}_{x}\\\
\end{array}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{2E}M^{2}\left[\begin{array}[]{c}{\nu}_{e}\\\
{\nu}_{x}\\\ \end{array}\right]$ (13) $\displaystyle=$
$\displaystyle\frac{1}{2E}\left[U\left[\begin{array}[]{cc}m_{1}^{2}&0\\\
0&m_{2}^{2}\\\ \end{array}\right]U^{\dagger}+\left[\begin{array}[]{cc}A&0\\\
0&0\\\ \end{array}\right]\right]\left[\begin{array}[]{c}{\nu}_{e}\\\
{\nu}_{x}\\\ \end{array}\right]$ (20) $\displaystyle=$
$\displaystyle\frac{1}{4E}\left[(\Sigma+A)+\left[\begin{array}[]{cc}{A-{\Delta}C_{2\theta}}&{\Delta}S_{2\theta}\\\
{\Delta}S_{2\theta}&{-A+{\Delta}C_{2\theta}}\\\
\end{array}\right]\right]\left[\begin{array}[]{c}{\nu}_{e}\\\ {\nu}_{x}\\\
\end{array}\right],$ (25)

where $\Sigma=m_{2}^{2}+m_{1}^{2}$, $\Delta=m_{2}^{2}-m_{1}^{2}$, the mixing
angle $\theta$ refers to $\theta_{13}$ and $C_{2\theta}=\cos{2\theta}$,
$S_{2\theta}=\sin{2\theta}$. In recent years, the new data about $\theta_{13}$
have been provided by some neutrino experiments Ardellier Abe1 -An1 Abe2 Ahn .
In the following, we will use $\theta_{13}\simeq 8.8^{\circ}$ provided by the
Daya Bay Collaboration An1 .

By solving Eq. (4), we can get how $\nu_{e}$ and $\nu_{x}$ evolute with
distance, then the survival probability and conversion probability are given
as follows:

$\displaystyle
P_{s}=P({\nu_{e}}\rightarrow{\nu_{e}})={\left|\left\langle{{\nu}_{e}}|{{\nu}_{e}}\right\rangle\right|}^{2},$
$\displaystyle
P_{c}=P({\nu_{x}}\rightarrow{\nu_{e}})={\left|\left\langle{{\nu}_{e}}|{{\nu}_{x}}\right\rangle\right|}^{2}.$
(26)

##  III The density profiles of different SN models

In the previous study, we used the simplified density profile to represent the
real state of density without considering the mass of the supernova Guo -Xu .
In this section, based on simulation data we can provide the density profiles
at different times in three supernova models. These three supernova models are
8.8, 10.8 and 18 times the mass of the sun, respectively, so they will be
called as the 8.8$M_{\odot}$ model, the 10.8$M_{\odot}$ model and the
18$M_{\odot}$ model respectively for simplicity in the following.

###  III.1 The 8.8$M_{\odot}$ model

In Fig. 1 it is shown how the density changes with the radius within 6 seconds
after the supernova explosion in the 8.8$M_{\odot}$ model.

It can be seen that on each density curve there is an obvious decrease at
about $10^{1.5}$ km of radius, and this decline gets more and more notable
with the increase of time. For example, at the time of 0.2 s, the density
declines about $10^{5}$ g/cm3 from 10 km to $10^{2}$ km, while when $t$ = 6 s,
the density declines about $10^{7}$ g/cm3 from 10 km to $10^{1.5}$ km. After
the sharp decrease, the density continues to decline, but relatively slowly.
In this slow-decline stage, from $10^{2}$ km - $10^{4}$ km, the density goes
down about $10^{7}$ g/cm3. Except for the curve of 1 s, at the end of other
curves, there is an abnormal increase of density. Compared with the rise with
some twists and turns on the curves of 2 s $\sim$ 6 s, at the end of 0.2 s and
0.5 s curves, the densities ascend slight increases. Yet on the curve of 1 s,
there is no obvious increase at all.

Overall, in the 8.8$M_{\odot}$ model, no matter what time it is, the density
decreases with the increase of radius if we do not take the ends of the curves
into account. It is worth noting that the so-called shock wave caused by the
changing density does not appear in the 8.8$M_{\odot}$ model at all. This is
remarkably different from the simplified density profile we used previously
Fogli2 .

###  III.2 The 10.8$M_{\odot}$ model

In Fig. 2, there are 11 curves to demonstrate how the density varies with the
radius at different times in the 10.8$M_{\odot}$ model.

Similar to Fig. 1, there is an abrupt density decrease at about $10^{1.5}$ km
of radius on each curve in Fig. 2. Then the density continues to decline
slowly in comparison with the previous sharp decrease. During this phase, it
can be observed that there is a density ”bulge” on each curve, which is the
forward shock wave we studied before. At 0.5 s, it seems that this shock wave
is a bit fuzzy, but after 1 s it becomes more and more evident and moves
forward with time. On the curves of $t\geqslant$ 3 s, it can be seen that
there is a density ”concave” located behind the forward shock wave. That is
the reverse shock wave, resulting from the supernova matter squeezing with the
end of the shock wave. Just like the forward shock wave, the reverse shock
wave also moves forward with time. On the curves of $t\geqslant$ 6 s, a small
density ”bump” appears on the forward shock wave, indicating that the changes
of matter within the forward shock wave are very complicated. This phenomenon
has never be mentioned before Fogli2 Schirato Huang .

Furthermore, it is found that all the curves in Fig. 2 converge to a line at
the end of themselves. Through the curve fitting method, it is found that this
converging line can be represented by $\rho=r^{-2.7}$ approximately, in close
proximity to $\rho=r^{-3}$, which is known as the expression of density
distribution in a stable supernova as mentioned in Refs. Schirato Kuo .

###  III.3 The 18$M_{\odot}$ model

Same as Fig. 2, Fig. 3 shows the density profiles at 11 different points of
time in the 18$M_{\odot}$ model. Just like the two former models, the density
in the 18$M_{\odot}$ model also decreases rapidly at about $10^{1.5}$ km
judged from the appearance of each curve and then the density declines
gradually. On the curve of $t$ = 0.5 s, it seems that the forward shock wave
is forming but not complete yet. While on other curves, it can be seen clearly
that there is a forward shock wave moving with time. Also at the end of all
curves, the density profiles of different times converge into a line. Compared
with Fig. 2, in which there is a density ”concave” on the curves of
$t\geqslant$ 3 s, it is found that the reverse shock wave is not formed
throughout in the 18$M_{\odot}$ model.

The above three figures of density profiles have something in common. That is,
at about $10^{1.5}$ km of radius, there is an obvious density decline which
gets more and more abrupt with time. It might be assumed that the core of the
supernova, which is incredibly dense, is collapsing more and more quickly with
time.

What is more, these three figures indicate that the mass of the supernova is
very closely related to the shape of shock wave. If the mass is not big
enough, taking the 8.8$M_{\odot}$ model as an example, there would be no shock
wave (neither the forward shock wave nor the reverse shock wave) appearing in
the supernova. However, if the supernova mass is too massive, for instance,
the 18$M_{\odot}$ model, the reverse shock wave would not arise in it. In
other words, only when the mass of the supernova is neither too small nor too
big, and within a certain range (from the three models, this range should be
$8.8M_{\odot}<$ the mass of SN $<18M_{\odot}$), both the forward shock wave
and the reverse shock wave will be generated in it. Moreover, the forward
shock wave and the reverse shock wave do not appear at the same time. The
forward shock wave is formed within one second after the supernova explosion,
while the reverse shock wave can been seen clearly after about 3 second of
explosion.

##  IV The variations of $P_{s}$ and $P_{c}$ in different SN models

In the previous section, the density profiles in three supernova models are
displayed graphically and we have already known that in the 8.8$M_{\odot}$
model shock wave never exist and in the 18$M_{\odot}$ model only the forward
shock wave can be generated, while both the forward and the reverse shock wave
are formed in the 10.8$M_{\odot}$ model during the explosion. No matter
whether the shock wave exist and no matter how its appearance is, as long as
we know the data of density profiles in supernova models, with the theory of
neutrino oscillations in matter as formulated in Section 2, we can calculate
the survival probability of $\nu_{e}\rightarrow\nu_{e}$, as well as the
conversion probability of $\nu_{x}\rightarrow\nu_{e}$. Given the discrete
density values at some points of radius, by substituting the density values
into Eq. (4) and through neutrinos’ point-by-point evolution, then we can
obtain the needed results.

###  IV.1 The 8.8$M_{\odot}$ model

In this subsection, based on the theory about the neutrino propagation in
matter and applying the density profiles in Section 3, we will show in detail
how $P_{s}$ and $P_{c}$ change in the 8.8$M_{\odot}$ supernova models,
respectively.

In Fig. 4, it is shown how the survival probability of $\nu_{e}$, $P_{s}$,
changes over time when the neutrino energy takes four different values.
Because of the time interval between the density profiles at different time
points, $P_{s}$ changes in a broken-line manner. As we can see from Fig. 4,
for $E$ = 11, 16 and 25 MeV, the value of $P_{s}$ has a significant reduction
from 0.2 s to 0.5 s. After 1 second, $P_{s}$ is changing all the time, but its
changes are very small, and the value of $P_{s}$ is basically no more than
0.05. While for the line of $E$ = 40 MeV, $P_{s}$ has small changes before 1 s
without dramatic decline, then with the increase of time its value twists and
turns up to 0.2 at 6 s.

Fig. 5 details how $P_{s}$ varies with time and the neutrino energy in a
three-dimensional graph. It looks that $P_{s}$ changes continuously because
the curves were processed by ”ployfit” under the assumption that there is no
sudden change between different time points. This method is also used in other
three-dimensional graphs in the following passage. It is found that at the
early stage of explosion ($t\leqslant$ 1 s) and when the neutrino energy is
relatively small ($E\leqslant$ 10 MeV), the value of $P_{s}$ is relatively
large but declines abruptly in a short time. These obvious changes form some
sharp protuberances at the corner of the figure. Then with the increase of
time and the neutrino energy, $P_{s}$ goes into the next stage, in which the
variations of $P_{s}$ are smoother compared with before. Although there are
many small fluctuations, the value of $P_{s}$ always keeps below 0.1 basically
in this stage. But as the time is later than 4 s and at high neutrino energy
($E\geqslant$ 40 MeV), some large fluctuations appear, which means $P_{s}$ is
greatly affected by time and the neutrino energy during this period. In
parallel with the increase of time and the neutrino energy, $P_{s}$ is rising
with fluctuations in the area of $t\geqslant$ 5 s and $E\geqslant$ 60 MeV, and
its value is even close to 0.5 ultimately.

At the same time, we also calculate the conversion probability of
$\nu_{x}\rightarrow\nu_{e}$, $P_{c}$, in the 8.8$M_{\odot}$ model. Fig. 6
shows the time-dependent changes of $P_{c}$, where the neutrino energy is
taken the same values as in Fig. 4. Different from $P_{s}$, the value of
$P_{c}$ is relatively large. When $t$ = 0.2 s, for $E$ = 11, 25 and 40 MeV,
$P_{c}$ is greater than 0.9; for $E$ = 16 MeV, $P_{c}$ is greater than 0.75.
With time increasing from 0.2 s to 0.5 s, the value of $P_{c}$ rises up to
about 0.97 for all the lines of different neutrino energy. In generally, for
$E$ = 11, 16, 25 MeV, $P_{c}$ always keeps above 0.95 in spite of some small
changes after 0.5 s. While for $E$ = 40 MeV, the value of $P_{c}$ shows a
downward trend roughly over the same time period, and decreases to about 0.82
at 6 s.

In Fig. 7, it can be seen how the $P_{c}$ varies with time and the neutrino
energy. At the beginning of explosion ($t\leqslant$ 1 s) and at low neutrino
energy ($E\leqslant$ 20 MeV), $P_{c}$ has abrupt changes so there are some
”upside-down peaks” in the corner of this figure. Then the value of $P_{c}$
goes up and down, and maintains between 0.8 and 1 basically. As time goes on
and the neutrino energy gets bigger, in the area of $t\geqslant$ 4 s and
$E\geqslant$ 50 MeV, the extent of variations become more and more remarkable.

By comparing Fig. 7 with Fig. 5, it seems that if $P_{s}$ can be rotated 180
degrees with respect to the XY plane, we would obtain the figure of $P_{c}$.
If so, the sum of $P_{s}$ and $P_{c}$ would equal 1, then it means the number
of $\nu_{e}$ is unchangeable even neutrino oscillations exist in the
supernova. If the result of $P_{s}$ added to $P_{c}$ is greater than 1, it
means the number of $\nu_{e}$ is more than that when they are emitted.
Otherwise, the result of less than 1 implies that the number of $\nu_{e}$ is
less than the original number. In order to find whether the number of
$\nu_{e}$ loses or increases, the results of $P_{s}$ added to $P_{c}$ as the
function of the neutrino energy are shown in Fig. 8.

Because the result of 0.2 s is quite a bit different from those of other
points of time, the results are shown in two figures. In the left figure, it
is found that when $t$ = 0.2 s, the number of $\nu_{e}$ is severely changeable
with the neutrino energy, which means at the beginning of explosion, the
conversion between $\nu_{e}$ and $\nu_{x}$ is very drastic. In the right
figure, we can see the results are always changing but trend toward stability
compared with the case of 0.2 s and keep within the range of 0.93 - 1.07. That
is, the losing and increasing ratio of $\nu_{e}$ are both no more than 7$\%$
during 0.5 s - 6 s of explosion.

###  IV.2 The 10.8$M_{\odot}$ model

In Fig. 9, there are four typical examples of $P_{s}$ at different neutrino
energies changing with time in the 10.8$M_{\odot}$ model. It is found that the
value of $P_{s}$ is less than 0.2 at 0.5 s for $E$ = 11, 16, 25, 40 MeV. When
$t$ = 1 s, $P_{s}$ of four different neutrino energies all drop to about 0.02
and keep the value till $8s$. As for $E$ =11, 16, 25 MeV, $P_{s}$ maintains
its value until 10 s; but for $E$ = 40 MeV, at 9 s $P_{s}$ increases to about
0.08 and then falls back to approach the previous value before 8 s. Compared
with Fig. 4, it seems that the time-varying fluctuations of $P_{s}$ in the
10.8$M_{\odot}$ model are much less than those in the 8.8$M_{\odot}$ model.

The variations of $P_{s}$ with time and the neutrino energy in the
10.8$M_{\odot}$ model is shown by a three-dimensional graph in Fig. 10. At 0.5
s, the value of $P_{s}$ changes significantly with the increase of the
neutrino energy, so there are many fluctuations along the axis of the neutrino
energy. From $t$ = 0.5 s to $t$ = 1 s, $P_{s}$ has an abrupt decrease, which
is in accord with the case shown in Fig. 9. With time going on, despite the
neutrino energy is different, $P_{s}$ goes into a relatively stable state and
its value is always around zero. But in the region of $t\geqslant$ 8 s and
$E\geqslant$ 40 MeV, some obvious fluctuations of $P_{s}$ revive again,
however, its value never exceeds 0.2.

Overall, there is an obvious difference between the 10.8$M_{\odot}$ model and
the 8.8$M_{\odot}$ model about $P_{s}$. That is, $P_{s}$ is more stable in the
10.8$M_{\odot}$ model than in the 8.8$M_{\odot}$ model. At the same time, two
models also have some similarities: at the beginning of explosion
($t\leqslant$ 1 s), $P_{s}$ has abrupt changes in a short time and fluctuates
significantly with the neutrino energy; in the area of large neutrino energy
and during the later time of explosion, fluctuations of $P_{s}$ become obvious
and frequent.

In Fig. 11, it can be seen how $P_{c}$ changes with time when the neutrino
energy takes four different values in the 10.8$M_{\odot}$ model. Obviously,
for all the lines the value of $P_{c}$ increases from more than 0.82 of 0.5 s
to about 0.98 of 1 s. Afterwards, for $E$ = 11, 16 MeV, $P_{c}$ keeps
invariable basically from 2 s to 9 s and decreases slightly at 10 s. On the
lines of $E$ = 25, 40 MeV, there are relatively small drops at 3 s and
relatively large drops at 9 s respectively, yet without significant changes at
other points of time.

The variations of $P_{c}$ with time and the neutrino energy in the
10.8$M_{\odot}$ model is shown in Fig. 12. Similar to $P_{s}$ in Fig. 10, at
0.5 s, $P_{c}$ has many fluctuations along the axis of the neutrino energy and
a sudden change from 0.5 s to 1 s. When $t\geqslant$ 8 s and $E\geqslant$ 40
MeV, the fluctuations of $P_{c}$ exist again, and its value varies between 0.8
and 1 approximately. While in other areas, $P_{c}$ maintains the basic
stability. Generally, the figure of $P_{c}$ in the 10.8$M_{\odot}$ model is
smoother by comparison with $P_{c}$ of the 8.8$M_{\odot}$ model shown in Fig.
7 .

Likewise, we also plot the results of $P_{s}$ added to $P_{c}$ to find how the
number of $\nu_{e}$ changes in the 10.8$M_{\odot}$ model. Same as Fig. 8,
there are two figures in Fig. 13. In the left figure, it can be seen that at
0.5 s, the results hover drastically around 1 with the increase of the
neutrino energy, which means at this point of time, the conversion of
$\nu_{e}$ and $\nu_{x}$ is very dramatic and greatly affected by the neutrino
energy. By contrast, the results of other points of time are relatively
stable. The right figure indicate that from 1 s to 10 s of the explosion, the
sum of $P_{s}$ and $P_{c}$ fluctuates in the range of 0.96 and 1.01. At 1 s,
the results exceed 1 occasionally, while from 2 s to 10 s, the results are
almost under the lever of 1, which illustrates that the number of $\nu_{e}$ is
more likely to lose during this period. The losing rate and the increasing
rate of $\nu_{e}$ is less than 4$\%$ and 1$\%$, respectively.

###  IV.3 The 18$M_{\odot}$ model

In Fig. 14, it is illustrated how $P_{s}$ changes as the function of time at
four different neutrino energies in the 18$M_{\odot}$ model. In this figure,
we can see that at 1 s, for the four neutrino energies, $P_{s}$ drops below
0.05 from its original value, and then it remains so until the end of time.

In Fig. 15, it is found that $P_{s}$ has many prominent fluctuations at 0.5 s
with the change of the neutrino energy. After 1 s, no change can be observed
so the value of $P_{s}$ almost becomes a plane in the three-dimensional graph.
This phenomenon corresponds to the straight lines from 1 s to 10 s in Fig. 13.
The value of $P_{s}$ is very close to 0, meaning that $\nu_{e}$ almost
completely convert into $\nu_{x}$ during this period.

In Fig. 16, for the four lines of different neutrino energies, the value of
$P_{c}$ in the 18$M_{\odot}$ model has an obvious rise at 1 s and then
basically remains constant although the time is changing. The only noticeable
one is a slight decline at 1 s for $E$ = 11 and 16 MeV.

In Fig. 17, it can be seen that $P_{c}$ changes with the neutrino energy very
significantly at 0.5 s. Just like Fig. 7 and Fig. 12, there are many
fluctuations on the edge of the figure. From 0.5 s to 1 s, the value of
$P_{c}$ goes up obviously, which can be seen more clearly in Fig. 16. After
the marked change, $P_{c}$ is nearly invariable and remains broadly flat over
the rest of the region consisting of time and the neutrino energy. And the
value of $P_{c}$ always near 1 suggests that almost all $\nu_{x}$ convert into
$\nu_{e}$ after 1 second of explosion.

In the left picture of Fig. 18, the line of 0.5 s varies with the neutrino
energy remarkably. But on other lines, it is difficult to observe significant
changes. In the right figure, it is shown that the variations of results at
low energies ($E\leqslant$ 30 MeV) are a little remarkable than those at high
energies during the time of 1 s - 10 s. But the change interval is always
between 0.975 and 1.01. In other words, the permissible losing rate of
$\nu_{e}$ is no more than 2.5$\%$ and the increasing rate is no more that
1$\%$ during this time.

##  V Summery and Discussions

By comparing the variations of $P_{s}$ in three models, it is observed that at
the start of explosion ($t<$ 1 s), $P_{s}$ changes abruptly in a short time
and changes significantly with the neutrino energy. Later on, it remains
relatively stable. But with the increase of time and the neutrino energy,
fluctuations of $P_{s}$ start to get more and more prominent in the
8.8$M_{\odot}$ model and the 10.8$M_{\odot}$ model. Meanwhile, the intensity
and the scale of fluctuations in the 10.8$M_{\odot}$ model are lower and
smaller than those in the 8.8$M_{\odot}$ model. Nevertheless, this ”restart-
fluctuating” phenomenon does not appear in the 18$M_{\odot}$ model, in
consequence the previous stable state is maintained until the end. In brief,
$P_{s}$ is getting more stable with the increase of supernova masses in the
models.

Similar to the variations of $P_{s}$, $P_{c}$ in three models also have an
abrupt change before 1 s and can be affected obviously by the neutrino energy
at the initial stage of explosion. Followed by a relatively stable phase,
$P_{c}$ restarts to fluctuate over time and with the increase of the neutrino
energy. Same as the case of $P_{s}$, with the increase of supernova mass in
the model, the relative stable area of $P_{c}$ is getting bigger and bigger.

By plotting the results of $P_{s}$ added to $P_{c}$ in three models, we can
analyze how the number of $\nu_{e}$ changes with the neutrino energy at
different times. It turns out that the number of $\nu_{e}$ is extremely
unstable at the beginning of explosion ($t<$ 1 s), reflecting the dramatic
conversions between $\nu_{e}$ and $\nu_{x}$. During the next period of
explosion (1 s - 6 s for the 8.8$M_{\odot}$ model; 1 s - 10 s for the
10.8$M_{\odot}$ and the 18$M_{\odot}$ model), the number of $\nu_{e}$ has some
small changes with time and the neutrino energy. Generally, if we could detect
$\nu_{e}$ after one second of explosion outside the supernova, compared with
the moment when $\nu_{e}$ is just emitted, the number of $\nu_{e}$ might be
increasing or reduced. Its change rate is no more than 7$\%$, 4$\%$ and
2.5$\%$ in the 8.8$M_{\odot}$, 10.8$M_{\odot}$ and 18$M_{\odot}$ models,
respectively. It seems that the number of $\nu_{e}$ is more stable in the more
massive supernova, although the conversions between $\nu_{e}$ and $\nu_{x}$
are very drastic.

In summary, the supernova mass in the model determines the density profiles in
the supernova. Consequently, it can affect the interactions between neutrinos
and matter, and then make changes to $P_{s}$, $P_{c}$ and the number of
neutrinos that can be detected. So it needs to be stressed that the supernova
mass is an important factor when properties of supernova neutrinos are
studied.

There is no doubt that if we want to obtain the number of neutrinos that can
be detected on the earth, discussing the shock effects is merely not enough.
In Refs. Guo Huang , we have already studied on the earth effects
comprehensively and the relevant researches on the collective effects are
ongoing currently Dasgupta -Hu . So far, we have already known that the
collective effects are very important on the conversions between different
neutrinos and affect other aspects about neutrinos, such as the mass hierarchy
and luminosity Chakraborty1 Fogli02 . With the progress of study on the
impacts of shock wave effects and collective effects, it is possible to obtain
the results in good agreement with experiments.

##  VI Acknowledgments

We are very grateful to Tobias Fisher for providing the simulation data about
the density profiles of supernovae. This work was supported in part by
National Science Foundation of China (Project Numbers 11175020, 11275025).

## References

  * (1) W. D. Arnett, Astrophys. J. A 319 136 (1987).
  * (2) R. M. Bionta, G. Blewitt, C. B. Bratton, D. Casper and A. Ciocio, et al., Phys. Rev. Lett. 58 1494 (1987).
  * (3) H. A. Bethe, Rev. Mod. Phys. 62 801 (1990); G. E. Brown, H. A. Bethe, and G. Baym, Nucl. Phys. A 375 481 (1982).
  * (4) F. Ardellier, I. Barabanov, et al., [arXiv:hep-ex/0606025].
  * (5) X.-H. Guo, et al. [Daya-Bay Collaboration], [arXiv:hep-ex/0701029].
  * (6) K. Abe, et al., [T2K Collaboration], Phys. Rev. Lett. 107 041801 (2011).
  * (7) P. Adamson, et al., [MINOS Collaboration], Phys. Rev. Lett. 107181802 (2011) .
  * (8) F. P. An, et al., [Daya-Bay Collaboration] Phys. Rev. Lett. 108 171803 (2012).
  * (9) F. P. An, et al., [Daya-Bay Collaboration] Nucl. Instrum. Meth. A 685 1 (2012).
  * (10) F. P. An, et al., [Daya-Bay Collaboration] Phys. Rev. Lett. 112 061801 (2014).
  * (11) Y. Abe, et al., [DOUBLE-CHOOZ Collaboration] Phys. Rev. Lett. 108 131801 (2012).
  * (12) J. K. Ahn, et al., [RENO Collaboration], Phys. Rev. Lett. 108, 191802 (2012)
  * (13) S. W. Bruenn, et al., [arXiv:1002.4909 [astro-ph.SR]].
  * (14) H. T. Janka, Ann. Rev. Nucl. Part. Sci. 62, 407 (2012) [arXiv:1206.2503 [astro-ph.SR]].
  * (15) S. Akiyama, J. C. Wheeler, D. L. Meier, I. Lichtenstadt, Astrophys. J. 584 954 (2003).
  * (16) K. Kotake, K. Sato, and K. Takahashi, Rept. Prog. Phys. 69 971 (2006).
  * (17) L. Hudepohl, et al., Phys. Rev. Lett. 104 251101 (2010).
  * (18) S. Chakraborty, et al., Phys. Rev. D 84 025002 (2011); Phys. Rev. Lett. 107 151101 (2011).
  * (19) D. N. Spergel, et al., Science 237 1471 (1987).
  * (20) C. Giunti and C. W. Kim, Oxford, UK: Univ. Pr. (2007) 710p.
  * (21) G. L. Fogli, et al., J. Cosmol. Astropart. Phys. 0504 002 (2005).
  * (22) K. Takahashi, K. Sato, H. E. Dalhed and J. R. Wilson, Astropart. Phys. 20, 189 (2003) [astro-ph/0212195].
  * (23) G. L. Fogli and E. Lisi, A. Mirizzi and D. Montanino, Phys. Rev. D 68 033005 (2003).
  * (24) R. C. Schirato and G. M. Fuller, Preprint [astro-ph/0205390] (2002).
  * (25) K. Takahashi, et al., Astropart. Phys. 20 189 (2003).
  * (26) R. Tomas, et al., Astropart. Phys. 0409 015 (2004).
  * (27) L. Wolfenstein, Phys. Rev. D 17 2369 (1978); 20 2634 (1979).
  * (28) T. K. Kuo and J. Pantaleone, Rev. Mod. Phys. 61 937 (1989).
  * (29) S. P. Mikheyev and A. Y. Smirnov, Nucl. Phys. B 42 913 (1985).
  * (30) M. Kachelriess, A. Strumia, R. Tomas, J. W. F. Valle, Phys. Rev. D 65 073016 (2002).
  * (31) X.-H. Guo, M.-Y. Huang, and B.-L. Young, Phys. Rev. D 79 113007 (2009).
  * (32) M.-Y. Huang, X.-H. Guo, and B.-L. Young, Phys. Rev. D 82 033011 (2010).
  * (33) XU Jing, HUANG Ming-Yang, HU Li-Jun, Guo Xin-Heng and YOUNG Bing-Lin, Commun. Theor. Phys. 61 226 (2014).
  * (34) J. Kneller, C. Volpe, Phys. Rev. D 82 123004 (2010).
  * (35) T. Lund and J. P. Kneller, Phys. Rev. D 88 02308 (2013).
  * (36) S. Galais, et al., Phys. Rev. D 81 053002 (2010).
  * (37) K. M. Patton, J. P. Kneller and G. C. McLaughlin, Phys. Rev. D 89 073022 (2014). [arXiv:1310.5643 [hep-ph]].
  * (38) A. Wongwathanarat, E. Mueller and H.-T. Janka, [arXiv:1409.5431 [astro-ph.HE]].
  * (39) S. M. Couch and C. D. Ott, [arXiv:1408.1399 [astro-ph.HE]].
  * (40) B. Dasgupta and A. Dighe, Phys. Rev. D 77 113002 (2008).
  * (41) S. Chakraborty and S. Choubey, et al., J. Cosmol. Astropart. Phys. 1006 007 (2010).
  * (42) H. Y. Duan and J. P. Kneller, J. Phys. G 36 113201 (2009).
  * (43) H. Y. Duan, G. M. Fuller and Y. Z. Qian, Phys. Rev. D 74 123004 (2006); 76 085013 (2007); 75 125005 (2007); Phys. Rev. Lett. 99 241802 (2007).
  * (44) S. Hannestad, G. G. Raffelt, G. Sigl and Y. Y. Y. Wong, Phys. Rev. D 74 105010 (2006); 76 029901(E) (2007).
  * (45) G. G. Raffelt and A. Y. Smirnov, Phys. Rev. D 76 081301(R) (2007); 76 125008 (2007).
  * (46) G. L. Fogli and E. Lisi, et al., J. Cosmol. Astropart. Phys. 0712 010 (2007); 0910 002 (2009).
  * (47) Li-Jun Hu, Jing Xu, Xin-Heng Guo and Bing-Lin Young, in preparation.

Fig. 1

Figure 1: Density profiles of different times in the 8.8$M_{\odot}$ model. The
dark solid curve corresponds to 0.2 s. The blue solid, dashed, dotted and dot-
dashed curves correspond to 0.5 s, 1 s, 2 s and 3 s, respectively. The red
solid, dashed and dotted curves correspond to 4 s, 5 s and 6 s, respectively.

Fig. 2

Figure 2: Density profiles of different times in the 10.8$M_{\odot}$ model.
The dark dot-dashed curve corresponds to 0.5 s. The blue solid, dashed, dotted
and dot-dashed curves correspond to 1 s, 2 s, 3 s and 4 s, respectively. The
red solid, dashed, dotted and dot-dashed curves correspond to 5 s, 6 s, 7 s
and 8 s, respectively. The dark solid and dashed curves correspond to 9 s and
10 s, respectively.

Fig. 3

Figure 3: Density profiles of different times in the 18$M_{\odot}$ model. The
dark dot-dashed curve corresponds to 0.5 s. The blue solid, dashed, dotted and
dot-dashed curves correspond to 1 s, 2 s, 3 s and 4 s, respectively. The red
solid, dashed, dotted and dot-dashed curves correspond to 5 s, 6 s, 7 s and 8
s, respectively. The dark solid and dashed curves correspond to 9 s and 10 s,
respectively.

Fig. 4

Figure 4: $P_{s}$ as a function of the time $t$ for four different neutrino
energies in the 8.8$M_{\odot}$ model. The solid, dashed, dotted and dot-dashed
curves correspond to $E$ = 11, 16, 25 and 40 MeV, respectively. The cross
points are obtained from our calculations.

Fig. 5

Figure 5: The variations of $P_{s}$ with time and the neutrino energy in the
8.8$M_{\odot}$ model. The scopes of independent variables are: 0.2 s $\leq
t\leq$ 6 s, 1 MeV $\leq E\leq$ 80 MeV, respectively.

Fig. 6

Figure 6: $P_{c}$ as a function of the time $t$ for four different neutrino
energies in the 8.8$M_{\odot}$ model. The solid, dashed, dotted and dot-dashed
curves correspond to $E$ = 11, 16, 25 and 40 MeV, respectively. The cross
points are obtained from our calculations.

Fig. 7

Figure 7: The variations of $P_{c}$ with time and the neutrino energy in the
8.8$M_{\odot}$ model. The scopes of independent variables are: 0.2 s $\leq
t\leq$ 6 s, 1 MeV $\leq E\leq$ 80 MeV, respectively.

Fig. 8

Figure 8: The results of $P_{s}$ added to $P_{c}$ as a function of the
neutrino energy at different times in the 8.8$M_{\odot}$ model. The blue
solid, dashed, dotted and dot-dashed curves correspond to 0.2 s, 0.5 s, 1 s
and 2 s, respectively. The red solid, dashed, dotted and dot-dashed curves
correspond to 3 s, 4 s, 5 s and 6 s, respectively. In the left figure, there
are 8 curves of all different times. In the right figure, there are 7 curves
except the curve of 0.2 s.

Fig. 9

Figure 9: $P_{s}$ as a function of the time $t$ for four different neutrino
energies in the 10.8$M_{\odot}$ model. The solid, dashed, dotted and dot-
dashed curves correspond to $E$ = 11, 16, 25 and 40 MeV, respectively. The
cross points are obtained from our calculations.

Fig. 10

Figure 10: The variations of $P_{s}$ with time and the neutrino energy in the
10.8$M_{\odot}$ model. The scopes of independent variables are: 0.5 s $\leq
t\leq$ 10 s, 1 MeV $\leq E\leq$ 80 MeV, respectively.

Fig. 11

Figure 11: $P_{c}$ as a function of the time $t$ for four different neutrino
energies in the 10.8$M_{\odot}$ model. The solid, dashed, dotted and dot-
dashed curves correspond to $E$ = 11, 16, 25 and 40 MeV, respectively. The
cross points are obtained from our calculations.

Fig. 12

Figure 12: The variations of $P_{c}$ with time and the neutrino energy in the
10.8$M_{\odot}$ model. The scopes of independent variables are: 0.5 s $\leq
t\leq$ 10 s, 1 MeV $\leq E\leq$ 80 MeV, respectively.

Fig. 13

Figure 13: The results of $P_{s}$ added to $P_{c}$ as a function of the
neutrino energy at different times in the 10.8$M_{\odot}$ model. The blue dot-
dashed curve correspond to 0.5 s. The dark solid, dashed, dotted and dot-
dashed curves correspond to 1 s, 2 s, 3 s and 4 s, respectively. The red
solid, dashed, dotted and dot-dashed curves correspond to 5 s, 6 s, 7 s and 8
s, respectively. The blue solid and dashed curves correspond to 9 s and 10 s,
respectively. In the left figure, there are 11 curves of all different times.
In the right figure, there are 10 curves except the curve of 0.5 s.

Fig. 14

Figure 14: $P_{s}$ as a function of the time $t$ for four different neutrino
energies in the 18$M_{\odot}$ model. The solid, dashed, dotted and dot-dashed
curves correspond to $E$ = 11, 16, 25 and 40 MeV, respectively. The cross
points are obtained from our calculations.

Fig. 15

Figure 15: The variations of $P_{s}$ with time and the neutrino energy in the
18$M_{\odot}$ model. The scopes of independent variables are: 0.5 s$\leq
t\leq$ 10 s, 1 MeV $\leq E\leq$ 80 MeV, respectively.

Fig. 16

Figure 16: $P_{c}$ as a function of the time $t$ for four different neutrino
energies in the 18$M_{\odot}$ model. The solid, dashed, dotted and dot-dashed
curves correspond to $E$ = 11, 16, 25 and 40 MeV, respectively. The cross
points are obtained from our calculations.

Fig. 17

Figure 17: The variations of $P_{c}$ with time and the neutrino energy in the
18$M_{\odot}$ model. The scopes of independent variables are: 0.5 s $\leq
t\leq$ 10 s, 1 MeV $\leq E\leq$ 80 MeV, respectively.

Fig. 18

Figure 18: The results of $P_{s}$ added to $P_{c}$ as a function of the
neutrino energy at different times in the 18$M_{\odot}$ model. The blue dot-
dashed curve correspond to 0.5 s. The dark solid, dashed, dotted and dot-
dashed curves correspond to 1 s, 2 s, 3 s and 4 s, respectively. The red
solid, dashed, dotted and dot-dashed curves correspond to 5 s, 6 s, 7 s and 8
s, respectively. The blue solid and dashed curves correspond to 9 s and 10 s,
respectively. In the left figure, there are 11 curves of all different times.
In the right figure, there are 10 curves except the curve of 0.5 s.