††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Ice front shaping by upward convective current

Ziqi Wang  Center for Combustion Energy, Key laboratory for Thermal Science
and Power Engineering of Ministry of Education, Department of Energy and Power
Engineering, Tsinghua University, Beijing 100084, China  Linfeng Jiang  Center
for Combustion Energy, Key laboratory for Thermal Science and Power
Engineering of Ministry of Education, Department of Energy and Power
Engineering, Tsinghua University, Beijing 100084, China  Yihong Du  Center for
Combustion Energy, Key laboratory for Thermal Science and Power Engineering of
Ministry of Education, Department of Energy and Power Engineering, Tsinghua
University, Beijing 100084, China  Chao Sun  Center for Combustion Energy, Key
laboratory for Thermal Science and Power Engineering of Ministry of Education,
Department of Energy and Power Engineering, Tsinghua University, Beijing
100084, China  Department of Engineering Mechanics, School of Aerospace
Engineering, Tsinghua University, Beijing 100084, China  Enrico Calzavarini
Univ. Lille, Unité de Mécanique de Lille - J. Boussinesq - UML - ULR 7512,
F-59000 Lille, France

###### Abstract

The extent and the morphology of ice forming in a differentially heated cavity
filled with water are studied by means of experiments and numerical
simulations. We show that the main mechanism responsible for the ice shaping
is the existence of a cold upward convective current in the system. Such a
current is ascribed to the peculiar equation of state of water, i.e., the non-
monotonous dependence of density with temperature. The precise form of the ice
front depends on several factors, first the temperature difference across the
cell which drives the convection, second the wall inclination with respect to
the vertical, both of which are here explored. We propose a boundary-layer
model and a buoyancy-intensity model which account for the main features of
the ice morphology.

Turbulent convective flows along with the ice formation process create
intriguing coupling behaviors, which have a widespread appearance in nature
and strong relevance in industrial applications meakin2010geological ;
Alboussiere2010Melting ; epstein1983complex ; worster1997convection ;
scagliarini2020modelling . In general, the orientation of the temperature
gradient and the gravity vector are not parallel. Their angle can play an
important role in determining the ice front morphology and system heat
transfer performances. Examples are the surficial icing of lakes and rivers,
floating ice bodies (e.g., icebergs), ice bodies (e.g., ice shelf) extending
outward from the land into waters, and solidification in energy storage
technology huppert1978melting ; russell1980melting ; rignot2013ice ;
kousha2017effect ; pogorelova2019moving ; hester2021aspect . For water, the
coupled physics among the phase change, turbulent convection, and the density
anomaly (water density reaches a maximum, $\rho_{c}$, at the density peak
temperature, $T_{\text{c}}$ ($\approx 4^{\circ}$C)) bring more challenges: the
gravitationally stable and unstable stratification coexist, which can strongly
affect flow structures malm1998bottom ; kowalewski1999freezing ;
Veronis1963Penetrative ; arid2012numerical ; large2014penetrative ;
corcione2015penetrative ; ayurzana2016phase ; 2020_coupled ;
lecoanet2015numerical ; toppaladoddi2018penetrative ; dietsche1985influence ;
hu2017lattice .

In recent years, many studies are devoted to exploring the interplay between
convective flows and thermal stratification/phase-transition under different
system inclination. For penetrative convection Veronis1963Penetrative , some
of these have found that the inclination can induce the breakdown of the fluid
stratification due to density anomaly inaba1984natural ; quintino2017heat ;
quintino2018optimal . Others have explored the coupling dynamics of phase-
change and turbulent convection using phase-change material
kamkari2014experimental ; kamkari2017numerical ; zeng2017numerical ;
madruga2021effect . A recent work Wang  _et al._ (2021), in a freezing-from-
above system, showed that the density anomaly induced stratification has major
effects on the flow structures and the resulting ice front speed and
equilibrium state. However, there is still a lack of explorations of the
physical mechanisms behind the whole rich ice front morphology. One may ask:
How does the ice front morphology change when the system is tilted? What are
the hydrodynamical mechanisms that account for the extent and the complex ice
front morphology?

In this work, by combining experiments, numerical simulations, and theoretical
modeling, we aim to systematically explore the freshwater solidification and
its coupling with turbulent convective flows to understand the complex
behaviors of the ice front morphology at varying system inclination angle,
$\beta$ (unit: deg). The experiments are conducted in a classical Rayleigh-
Bénard (RB) convection system sig94 ; bod00 ; ahlers2009heat ; lohse2010small
; chi12  (a fluid layer confined between a cold top plate, with temperature,
$T_{t}$, and a hot bottom plate, with temperature, $T_{b}$), with a quasi two-
dimensional rectangular shape (aspect ratios $L_{x}/H=1$ and $L_{z}/H=1/4$,
with $H=24~{}cm$ in experiments). The working fluid is deionized ultrapure
water (Prandtl number Pr $\approx 11$). The simulations are performed by means
of the Ch4-project code Calzavarini (2019), which adopts a Lattice-Boltzmann
algorithm for the description of fluid and temperature dynamics, and an
enthalpy method for the ice evolution (details see the Supplemental Material
wang_suppl ) Succi (2001); Huber  _et al._ (2008); Rabbanipour Esfahani  _et
al._ (2018); Chen  _et al._ (2017); Faden  _et al._ (2019); Wang  _et al._
(2021). Since the water thermal expansion coefficient inverses at
$T_{\text{c}}$ ($\approx 4^{\circ}$C), here we use the non-monotonous
relationship of density with temperature for water near $T_{\text{c}}$ Gebhart
and Mollendorf (1977), $\rho(T)=\rho_{c}(1-\alpha^{*}|T-T_{\text{c}}|^{q})$,
with $\rho_{c}=999.972~{}kg/m^{3}$ the maximum density $T_{c}\approx
4^{\circ}$C, $\alpha^{*}=9.30\times 10^{-6}(K^{-q})$, and $q=1.895$. Both in
experiments and simulations $T_{t}$ is fixed at $-10^{\circ}$C. In the
simulations we neglect the microscopic physics leading to kinetic
undercooling, Gibbs-Thomson effect and the anisotropic growth/melting Dash
_et al._ (2006).

Figure 1: Comparison of the ice morphology at the equilibrium state in VC from
experiments (a, b, c), simulations with (d, e, g) and without (f) considering
the density anomaly. The heating condition is: $T_{b}=6^{\circ}$C (a, d),
$10^{\circ}$C (b, e, f), $12^{\circ}$C (c, g). (d-g) Temperature fields with
color, 0∘C (black line) and 4∘C (red line) isotherms, and velocity vectors
(black arrows). The red and blue arrows in panels (d-g) show the rotating
convective rolls.

In this study, we monitor the local and the global ice thickness, denoted
respectively $h_{i}(x,t)$ and $h_{i}(t)$ and expressed in units of the cell
height $H$ ($h_{i}=1$ means full solidification). We consider that the
equilibrium is reached when the standard deviation of $h_{i}(t)$ over a time-
window of about 8 minutes is less than 0.5%. We first show the comparison of
the ice front morphology at the equilibrium state under different heating
conditions, $T_{b}$, between experiments and simulations, under vertical
convection (VC, solidification from left with $\beta=90$). Fig. 1(a, b, c) are
the experimental results with $T_{b}=6^{\circ}$C, $10^{\circ}$C and
$12^{\circ}$C, respectively (more results see the Supplemental Material
wang_suppl ). The corresponding simulation results with properly considering
the density anomaly are shown in Fig. 1(d, e, g), representing a good
agreement with the experimental measurements. This indicates that the
simulation indeed can capture the correct behavior of the system. The ice
front morphology displays a drastic change as $T_{b}$ increases. This is due
to the competition of two convective rolls originating from the density
anomaly, i.e., the one originating from the hot plumes detaching from the hot
plate (red arrow in Fig. 1(d, e, g)) and the other from the cold upward
convective current along the ice front (blue arrow in Fig. 1(d, e, g)). The
strength of the rolls can be adjusted by changing $T_{b}$. When
$T_{b}=6^{\circ}$C (Fig. 1(a)), the whole ice is shielded by the upward
convective current, and thus the ice front is flat but only tilted. As $T_{b}$
increases to $10^{\circ}$C, the ice front becomes highly uneven and is thicker
when protected by the upward convective current, thinner at the bottom, and
thinnest at the top where the impingement of hot plumes increases the local
heat transfer. Further, when $T_{b}$ is even higher (Fig. 1(c, g)), the
anticlockwise roll is much stronger, so it is able to intensively penetrate
the clockwise roll and finally affects the ice front, resulting in thinner
averaged ice thickness but with a similar shape to that of
$T_{b}=10^{\circ}$C.

In the simulations, the flow structures and ice front profile are highly
different when neglecting the density anomaly (i.e., the density is a linear
function of the temperature, see Fig. 1(f)): there is only one convective
roll, and the ice front is flat and thinner at the top and thicker at the
bottom. The clear distinction between the simulation with (Fig. 1(e)) and
without (Fig. 1(f)) considering the density anomaly indicates that this
property is crucial to properly describe the ice formation in presence of
natural convection.

Figure 2: Boundary layer model to explain the ice front morphology in the VC
case. (a) Sketch of the model: the ice front at $T_{0}=0^{\circ}$C (black
line), the $T_{m}$ isotherm which is the outer boundary of the thermal
boundary (red dashed line), the $T_{c}$ isotherm (green dashed line). The
angle between the tangential direction of the ice front and the x-direction is
$\gamma$; the thickness normal to the ice front (the green thick line in the
ice) is $\frac{h_{i}(x)}{cos\gamma}$ (with $cos\gamma=\frac{dS}{dx}$). (b-d)
Comparison of ice front morphology among experiments (shaded area),
simulations (line) and the model (dashed line) for different $T_{b}$.  Figure
3: (a) Nu as a function of $\beta$. Insets: sketches of the arrangement of
stably- and unstably-stratified layers of RB regime, $\beta<90$ (left) and
flipped RB, $\beta>90$(right). Error bars (estimated based on the temporal
time series of Nu in one simulation) are smaller than the symbols. The way to
measure the systematic error (i.e., ensemble average, qualitatively reflected
by the scatter of the data) is to perform multiple and independent
simulations, which is numerically expensive and therefore is not considered
here. (b-h) Instantaneous temperature field in ice, $T_{i}(x,y)$, and water,
$T_{w}(x,y)$, and isotherms (black thin lines) at the equilibrium state for
different angles $\beta=[0,180]$. The red dashed line and red thick line in
panels (b-h) are 0∘C and 4∘C isotherms, respectively. (i) Global ice thickness
$h_{i}$ as a function of $\beta$. The shaded area shows the spatial variation
of $h_{i}(x)$ in the simulations. The parameters are $T_{b}=10^{\circ}$C and
$T_{t}=-10^{\circ}$C for all cases.

To better understand the morphology of the ice front at its steady state, we
introduce a simple model based on the idea that a developing thermal boundary
layer forms along the ice-water interface. We note that when the ice has
ceased to grow there exists a local balance between the heat flux across the
ice and the one across the external boundary layer. The intensity of such
fluxes can be estimated by considering the thermal difference across the ice
and the boundary layer and the specific geometry of the problem, as shown in
the sketch Fig. 2(a)). We introduce a curvilinear coordinate, $S$, measuring
the length of the ice front starting from its boundary point at the $x=0$
position, and that is linked to the local ice thickness by the arc-length
formula $S(x)=\int_{0}^{x}\sqrt{1+(\tfrac{d(h_{i}(\xi))}{d\xi})^{2}}\,d\xi$.
We now express heat flux balance in the direction normal to the ice surface
as,

$k_{i}\frac{(T_{0}-T_{t})}{h_{i}(x)\frac{dS(x)}{dx}}=k_{w}\frac{(T_{m}-T_{0})}{\delta_{T}(S(x))},$
(1)

where $k_{w}$ and $k_{i}$ are the thermal conductivity of water and ice,
$(T_{0}-T_{t})$ is the temperature difference in the ice and
$h_{i}(x)/\tfrac{dS(x)}{dx}$ is the ice thickness in the direction normal to
the ice front. The estimation of the heat-flux in the water involves making
further assumptions. First, we consider that the heat transport in the
boundary layer in the direction normal to the ice front is purely conductive,
and this is justified by the fact that for $Pr\gg 1$ the thermal boundary
layer is nested inside the viscous one sun2008experimental . Second, we assume
that the thermal difference across the boundary layer of the fluid layer is
half of the one in the adjacent recirculation region (denoted with the blue
arrow in Fig. 1(d, e, g)), i.e. $(T_{m}-T_{0})$, with $T_{m}=(T_{0}+T_{c})/2$.
Third, the boundary layer thickness $\delta_{T}(S)$ is assumed to vary along
the ice front, as in a developing vertical thermal boundary layer, with a
dependence that we take to be $\delta_{T}(S)=C_{1}(S+C_{2})^{1/4}$
bejan2013convection ; white2006viscous ; shishkina2016momentum . The latter
expression has two dimensionless parameters
$C_{1}=c\cdot(g[1-\rho(T_{m})/\rho_{c}]/(\nu\kappa))^{1/4}$ with $c$ the
proportional constant ($c\approx 5~{}\text{m}^{3/4}$, $c$ has the unit m3/4 to
make $C_{1}$ dimensionless) bejan2013convection ; white2006viscous ;
shishkina2016momentum , and the offset $C_{2}$ because of non-zero boundary
layer thickness at $x=0$, with
$C_{2}=(h_{i}(0)[k_{w}(T_{m}-T_{0})]/[k_{i}(T_{0}-T_{t})]/C_{1})^{1/4}$, with
the boundary ice thickness $h_{i}(0)$ as an input from the simulation results.
Since the sidewall is adiabatic, $\frac{dh_{i}}{dx}|_{x=0}=0$ is an extra
known condition. With the above choices Eq.(1) becomes an integrodifferential
equation that can be solved numerically for the local ice thickness
$h_{i}(x)$. Fig.2(b-d) show the comparison of the model prediction with
experiments and simulations. A good qualitative agreement is reached in the
region where the upward convective current takes place. The disagreement in
the upper part of the ice front is expected as the boundary layer no longer
develops in that region due to the downward warmer convective current (denoted
with the red arrow in Fig. 1(d, e, g)). Although the exact spatial dependence
of $\delta_{T}$ along the ice interface is not known, the present model put
forward a robust physical mechanism for the shaping of the ice in the bottom
part of the convection cell.

Next, we perform systematical numerical simulations to explore how $\beta$
affects the extent and morphology of the ice front, with 0$\leq\beta\leq 180$
(with coordinate system attached to the cell). We limit our study to
$T_{b}=10^{\circ}$C, and the results are easy to extend to other situations.

First, we calculate the heat transfer rate, which when expressed
dimensionlessly is the global Nusselt number (Fig. 3(a)), Nu$=(\left\langle
u_{y}T\right\rangle-\kappa\partial_{y}\left\langle
T\right\rangle)/(\kappa\Delta T/H)$, where $\langle...\rangle$ represents an
average over time and the whole cell volume, $H$ is the system height. It is
noteworthy that there is convection both at $\beta=0$ (heating from below) and
$\beta=180$ (heating from above). Remarkably, as $\beta=180$, the fluid is
unstably stratified in the temperature interval between $T_{0}$ and $T_{c}$
which accounts for creating the convection. The latter feature is specific to
water and does not occur in other systems with working fluid’s density
increasing with temperature.

Figure 4: (a) Sketch of the buoyancy-intensity model to explain the position
of the maximum ice thickness. The buoyancy intensity of the clockwise and
anticlockwise convective rolls are $I_{1}(\beta)$ and $I_{2}(\beta)$,
respectively. $\Delta x_{0}$ locates the position of the maximum of the
$h_{i}(x)$, which is dimensionless normalized by $H$. $h_{i}(0)$, an input
parameter of the boundary layer model, is roughly constant and around $h_{i}$.
$h_{i}(0)$ vs. $\beta$ is reported in the Supplemental Material wang_suppl .
(b) $\Delta x_{0}$ as a function of $\beta$and its comparison with model
predictions (adjustable constant C=0.08). (c-m) Comparison of the ice front
morphology between simulations (shaded area) and the model (dashed line) under
different $\beta$ near $90$.

Fig. 3(b-h) show the temperature field with different $\beta$ from simulations
(more results see the Supplemental Material wang_suppl ). The inclination
results in different levels of thermal stratification, which induces huge
modifications to the ice front morphology. When $\beta=0$, there is a stably-
stratified layer (from $T_{0}$ to $T_{c}$) on top of the unstably-stratified
layer (from $T_{c}$ to $T_{b}$). As the cell is tilted with a small $\beta$,
the convection is strong enough to squeeze the stably-stratified layer to be
closely attached to the ice front, so the ice morphology is influenced by a
single-roll convective flow (see Fig. 3(b, c)). As $\beta$ increases towards
$90$, the inclination of the temperature gradient with respect to the gravity
is strong enough to break down the stratification and thus the stably-
stratified layer is set into motion in the form of a clockwise convective roll
(an upward cold water current) which competes with the initially anticlockwise
roll (downward warm water current). The shielding effect is prominent because
the ice thickness reaches a local maximum, and the hot plumes impacting on the
top part of the ice induces a local minimum of the ice thickness, based on
these two kinds of effects, the ice front presents an inflection point during
the transition from the thickest to the thinnest part, which has also been
reported in Fig. 1(b, c). The flow motion of the original stably-stratified
layer intensifies as $\beta$ increases for $\beta<90$. Beyond $\beta=90$, the
stratification configuration flips over. On the whole, the intensity of
convection is higher at $0<\beta<90$ than that of $90<\beta<180$. This is to
be connected to the different thermal difference across the respective
unstably stratified layers which are $(T_{c}-T_{0})\approx 4$K in the former
and $(T_{b}-T_{c})\approx 6$K in the latter case. This change in driving
strength also accounts for the observed results oF Nu (Fig. 3(a)).

To account for the influence of $\beta$ in a quantitative way, we now
calculate stationary global ice thickness ratio, $h_{i}=\langle
h_{i}(x,y)\rangle_{x,t}$, where $\langle...\rangle_{x,t}$ represents an
average over time and x-axis direction. Fig. 3(i) reports $h_{i}$, as a
function of $\beta$. $h_{i}$ shows an increasing trend as $\beta$ increases
because the heat transfers less efficiently for large $\beta$ (i.e., flipped
RB regime). The results from experiments (red circles) and simulations (green
squares) agree well with each other. It is noteworthy that the ice front is
highly variable in space as a result of different coupling with turbulent flow
structures. Here, the spatial fluctuations of the ice thickness are
represented by the local maximum and minimum of $h_{i}$, which is highlighted
by shaded area.

Coming back to cases when $\beta$ is around 90, we can observe a peculiar yet
robust form of ice front morphology (which can be observed in the range of
$\beta=40\sim 140$). In fact, the aforementioned boundary layer model can be
extended to the tilted system, by modifying $C_{1}$, with
$C_{1}=((g_{x}(1-\rho(T_{m})/\rho_{c})-g_{y}\mathcal{H}[\beta-90^{\circ}](1-\rho(T_{0})/\rho_{c}))/(\nu\kappa))^{1/4}$,
where $\mathcal{H}$ is Heaviside step function, $g_{x}=gsin\beta$,
$g_{y}=gcos\beta$. The first term in $C_{1}$ results from the inclination
effect, the second part originates from the inherent buoyancy contribution
induced by the density difference, and is present only for inclinations larger
than $90^{\circ}$. Fig. 4(c-m) show the comparison of the ice front morphology
between simulations (shaded area) and model prediction (dashed line). As
shown, the model can qualitatively capture behaviors of the ice front at the
inception of the thermal boundary. Note when $\beta<40$ and $\beta>140$, the
boundary layer attaching to the ice front is disturbed because of the plume
impacting under intensive interactions of the stably- and unstably-stratified
layers, so the model cannot be utilized. A second limitation of the model is
that the adopted expression for $\delta_{T}(S)$ is based on VC
bejan2013convection  but it does not involve possible dependencies on the cell
inclination. Nevertheless, it is remarkable that the model already performs
well in a wide range of $\beta$. Further studies are needed to improve or
quantitatively refine the model.

Another feature of the ice front is the position, $\Delta x_{0}$, where
$h_{i}(x)$ reaches maximum (symbols in Fig. 4(b)) and also implies the extent
of hot plume impacting effect. As discussed before, the local maximum of
$h_{i}$ originates from the competition of the buoyancy intensity between two
counterrotating rolls. The buoyancy intensity of the cold (clockwise) roll
(blue arrow in Fig. 4(a)), $I_{1}(\beta)$, and of the warm (anticlockwise)
roll (red arrow in Fig. 4(a)), $I_{2}(\beta)$, can be approximately evaluated
with

$\begin{split}&I_{1}(\beta)=g_{x}(1-\tfrac{\rho(T_{m})}{\rho_{c}})-g_{y}\mathcal{H}[\beta-90^{\circ}](1-\tfrac{\rho(T_{0})}{\rho_{c}});\\\
&I_{2}(\beta)=g_{x}(1-\tfrac{\rho(T_{m2})}{\rho_{c}})+g_{y}\mathcal{H}[90^{\circ}-\beta](1-\tfrac{\rho(T_{b})}{\rho_{c}}).\\\
\end{split}$ (2)

where the mean temperature $T_{m}=(T_{0}+T_{c})/2$ and
$T_{m2}=(T_{b}+T_{c})/2$. Eqn. (S5) holds when $\beta\neq 90$ (details see the
Supplemental Material wang_suppl ). The intensity ratio,
$(I_{1}(\beta)/I_{2}(\beta))$, captures the trend of $\Delta x_{0}$ as a
function of $\beta$, at least for $10<|\beta-90|<50$. A quantitative agreement
is obtained by an adjustment multiplicative factor $C\simeq 0.08$ (red dashed
line in Fig. 4(b)). This heuristic model provides further evidence that the
competition of two convective rolls accounts for the form of ice morphology.

To summarize, we found that the existence of a cold upward convective current,
due to the density anomaly of water, accounts for the ice shaping. We provide
physical understanding on the main features of the ice morphology in a wide
range of system inclinations. The present exploration offers deeper insight
into comprehending the liquid-solid interface morphology induced by the
coupling between phase-transition and natural convection with possible
applications in geophysical and climate sciences.

###### Acknowledgements.

This work was supported by Natural Science Foundation of China under grant nos
11988102, 91852202 and 11861131005.

## References

  * (1) P. Meakin, B. Jamtveit, Geological pattern formation by growth and dissolution in aqueous systems, Proc. R. Soc. A 466, 659 (2010).
  * (2) T. Alboussiere, R. Deguen, M. Melzani, Melting-induced stratification above the earth’s inner core due to convective translation, Nature 466, 744 (2010).
  * (3) M. Epstein, F. Cheung, Complex freezing-melting interfaces in fluid flow, Annu. Rev. Fluid Mech. 15, 293 (1983).
  * (4) M. G. Worster, Convection in mushy layers, Annu. Rev. Fluid Mech. 29, 91 (1997).
  * (5) A. Scagliarini, E. Calzavarini, D. Mansutti, F. Toschi, Modelling sea ice and melt ponds evolution: sensitivity to microscale heat transfer mechanisms, in: Mathematical Approach to Climate Change and its Impacts, (Springer, 2020) pp. 179–198.
  * (6) H. E. Huppert, J. S. Turner, On melting icebergs, Nature 271, 46 (1978).
  * (7) D. Russell-Head, The melting of free-drifting icebergs, Ann. Glaciol. 1, 119 (1980).
  * (8) E. Rignot, S. Jacobs, J. Mouginot, B. Scheuchl, Ice-shelf melting around antarctica, Science 341, 266 (2013).
  * (9) N. Kousha, M. Hosseini, M. Aligoodarz, R. Pakrouh, R. Bahrampoury, Effect of inclination angle on the performance of a shell and tube heat storage unit–an experimental study, Appl. Therm. Eng. 112, 1497 (2017).
  * (10) A. V. Pogorelova, V. L. Zemlyak, V. M. Kozin, Moving of a submarine under an ice cover in fluid of finite depth, J. Hydrodyn. 31, 562 (2019).
  * (11) E. W. Hester, C. D. McConnochie, C. Cenedese, L.-A. Couston, G. Vasil, Aspect ratio affects iceberg melting, Phys. Rev. Fluids 6, 023802 (2021).
  * (12) J. Malm, Bottom buoyancy layer in an ice-covered lake, Water Resour Res 34, 2981 (1998).
  * (13) T. A. Kowalewski, M. Rebow, Freezing of water in a differentially heated cubic cavity, Int. J. Comput. Fluid. D. 11, 193 (1999).
  * (14) G. Veronis, Penetrative convection, Astrophys. J. 137, 641 (1963).
  * (15) A. Arid, T. Kousksou, S. Jegadheeswaran, A. Jamil, Y. Zeraouli, Numerical simulation of ice melting near the density inversion point under periodic thermal boundary conditions, Fluid Dynamics & Materials Processing 8, 257 (2012).
  * (16) E. Large, C. Andereck, Penetrative Rayleigh-Bénard convection in water near its maximum density point, Phys. Fluids 26, 094101 (2014).
  * (17) M. Corcione, A. Quintino, Penetrative convection of water in cavities cooled from below, Comput. Fluids. 123, 1 (2015).
  * (18) B. Ayurzana, T. Hosoyamada, Phase change simulations of water near its density inversion point by lattice boltzmann method, in: Proceedings of the 23rd IAHR International Symposium on ice, (2016) pp. 1–8.
  * (19) P. Léard, B. Favier, P. Le Gal, M. Le Bars, Coupled convection and internal gravity waves excited in water around its density maximum at $4^{\circ}$C, Phys. Rev. Fluids 5, 024801 (2020).
  * (20) D. Lecoanet, M. Le Bars, K. J. Burns, G. M. Vasil, B. P. Brown, E. Quataert, J. S. Oishi, Numerical simulations of internal wave generation by convection in water, Phys. Rev. E 91, 063016 (2015).
  * (21) S. Toppaladoddi, J. S. Wettlaufer, Penetrative convection at high Rayleigh numbers, Phys. Rev. Fluids 3, 043501 (2018).
  * (22) C. Dietsche, U. Müller, Influence of bénard convection on solid–liquid interfaces, J. Fluid Mech. 161, 249 (1985).
  * (23) Y. Hu, D. Li, S. Shu, X. Niu, Lattice boltzmann simulation for three-dimensional natural convection with solid-liquid phase change, Int. J. Heat Mass Transf. 113, 1168 (2017).
  * (24) H. Inaba, T. Fukuda, Natural convection in an inclined square cavity in regions of density inversion of water, J. Fluid Mech. 142, 363 (1984).
  * (25) A. Quintino, E. Ricci, S. Grignaffini, M. Corcione, Heat transfer correlations for natural convection in inclined enclosures filled with water around its density-inversion point, Int. J. Therm Sci. 116, 310 (2017).
  * (26) A. Quintino, E. Ricci, S. Grignaffini, M. Corcione, Optimal inclination for maximum convection heat transfer in differentially-heated enclosures filled with water near $4^{\circ}$C, Heat Transf. Eng. 39, 499 (2018).
  * (27) B. Kamkari, H. Shokouhmand, F. Bruno, Experimental investigation of the effect of inclination angle on convection-driven melting of phase change material in a rectangular enclosure, Int. J. Heat Mass Transf. 72, 186 (2014).
  * (28) B. Kamkari, H. J. Amlashi, Numerical simulation and experimental verification of constrained melting of phase change material in inclined rectangular enclosures, Int. Commun. Heat Mass Transf. 88, 211 (2017).
  * (29) L. Zeng, J. Lu, Y. Li, W. Li, S. Liu, J. Zhu, Numerical study of the influences of geometry orientation on phase change material’s melting process, Adv. Mech. Eng. 9, 1 (2017).
  * (30) S. Madruga, J. Curbelo, Effect of the inclination angle on the transient melting dynamics and heat transfer of a phase change material, Phys. Fluids 33, 055110 (2021).
  * (31) Z. Wang, E. Calzavarini, C. Sun, F. Toschi, How the growth of ice depends on the fluid dynamics underneath, Proc. Natl. Acad. Sci. 118 (10), e2012870118 (2021).
  * (32) E. D. Siggia, High Rayleigh number convection, Annu. Rev. Fluid Mech. 26, 137 (1994).
  * (33) E. Bodenschatz, W. Pesch, G. Ahlers, Recent developments in Rayleigh-Bénard convection, Ann. Rev. Fluid Mech. 32, 709 (2000).
  * (34) G. Ahlers, S. Grossmann, D. Lohse, Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection, Rev. Mod. Phys. 81, 503 (2009).
  * (35) D. Lohse, K.-Q. Xia, Small-scale properties of turbulent Rayleigh-Bénard convection, Annu. Rev. Fluid Mech. 42, 335 (2010).
  * (36) F. Chillà, J. Schumacher, New perspectives in turbulent Rayleigh-Bénard convection, Eur. Phys. J. E 35, 58 (2012).
  * (37) See Supplemental Material at URL for the introduction of experimental and numerical methods, systematic investigation of various inclination angles, the boundary layer model and the buoyancy-intensity model.
  * (38) E. Calzavarini, Eulerian–lagrangian fluid dynamics platform: The ch4-project, Software Impacts 1,100002 (2019).
  * (39) S. Succi, The Lattice-Boltzmann equation: for fluid dynamics and beyond (Oxford university press, 2001).
  * (40) C. Huber, A. Parmigiani, B. Chopard, M. Manga, O. Bachmann, Lattice-boltzmann model for melting with natural convection, Int. J. Heat & Fluid Flow 29, 1469 (2008).
  * (41) B. Rabbanipour Esfahani, S. C. Hirata, S. Berti, E. Calzavarini, Basal melting driven by turbulent thermal convection, Phys. Rev. Fluids 3, 053501 (2018).
  * (42) S. Chen, Y. Y. Yan, and W. Gong, A simple lattice Boltzmann model for conjugate heat transfer research, Int. J. Heat and Mass Transf. 107, 862 (2017)
  * (43) M. Faden, A. Knig-Haagen, and D. Brggemann, An Optimum Enthalpy Approach for Melting and Solidification with Volume Change, Energies 12, 868 (2019).
  * (44) B. Gebhart, J. C. Mollendorf, A new density relation for pure and saline water, Deep Sea Research 24, 831 (1977).
  * (45) J. Dash, A. Rempel, J. Wettlaufer, The physics of premelted ice and its geophysical consequences, Rev. Mod. Phys. 78, 695 (2006).
  * (46) C. Sun, Y.-H. Cheung, K.-Q. Xia, Experimental studies of the viscous boundary layer properties in turbulent Rayleigh-Bénard convection, J. Fluid Mech. 605, 79 (2008).
  * (47) A. Bejan, Convection heat transfer (John wiley & sons, 2013).
  * (48) F. M. White, I. Corfield, Viscous fluid flow, Vol. 3 (McGraw-Hill New York, 2006).
  * (49) O. Shishkina, Momentum and heat transport scalings in laminar vertical convection, Phys. Rev. E 93, 051102 (2016).

Supplementary Information for:
Ice front shaping by upward convective current

##  Appendix A Experimental Methods

The experiments of freshwater solidification in a convective cell with
different inclination angles (here we use the coordinate system attached to
the cell) are conducted in a classical Rayleigh-Bénard (RB) convection system
(see Fig. S1). Fig. S1(a) shows the sketch of the experimental cell. The
experimental cell consists of a top cooling plate (with temperature $T_{t}$),
a heating bottom plate (with temperature $T_{b}$), and plexiglas sidewalls
with height H = 240 mm (length $L_{x}=240$ mm and width $L_{z}=60$ mm, i.e.,
aspect ratio $\Gamma=L_{x}/H=$1.0). The working fluid, which is confined
between the top and bottom plates, is deionized and ultrapure water, and the
water density and thermal expansion coefficient around the density peak
temperature, $T_{c}=$4∘C, are shown in Fig. S1(d) and (e). Before conducting
any experiments, water is boiled twice to be degassed. The temperature of both
the top and bottom plates are well controlled at constant temperatures
($T_{t}<0^{\circ}$C, and $T_{b}>0^{\circ}$C) by the circulating bath
(PolyScience PP15R-40), with the temperature fluctuations less than $\pm
0.2$K. Silicone O-rings are attached in between the top plate and the
sidewall, and also in between the sidewall and the bottom plate, to seal the
cell. In order to compensate the volume change during the solidification of
freshwater, an expansion vessel is connected to the experimental cell through
a rubber tube, which is open to the atmosphere so that the pressure of the
experimental cell remains unchanged. Six resistance thermistors (44000 series
thermistor element, see the sketch in Fig. S1(c)) are embedded into the top
and bottom plates, respectively (see the sketch in Fig. S1(a), the black dots
on the top and bottom plates show the placement of the thermistors). There are
two steps to avoid the heat exchange between the experimental cell and the
surrounding environment, 1) the experimental cell is wrapped in a sandwich
structure: insulation foam, aluminum plate (helping to support and acting as
the temperature measuring spot of the temperature sensor, PT 100, of the
Proportional-Integral-Derivative (PID) system), and insulation foam; 2) a PID
(Proportional-Integral-Derivative) controller is used to control the
temperature of the surrounding environment outside the experimental cell (see
Fig. S1(b)).

Figure S1: (a) Sketch of the experimental system for Rayleigh-Bénard
convection coupled with solidification of fresh water under different
inclination angles. (b)The PID (Proportional-Integral-Derivative) controller
and the temperature sensor used in the PID system. (c) The sketch of
resistance thermistor, 44000 series thermistor element, which is used to
measure the top and bottom plates temperature. (d) The nonmonotonic
relationship of the water density with temperature for cold water near $T_{c}$
from Ref. Gebhart and Mollendorf (1977). (e) the thermal expansion coefficient
of water as a function of the temperature.

During the experiments and numerical simulations, in order to focus on the
effect of inclination, we have limited the cold plate to a constant
undercooling temperature of $T_{t}=-10^{\circ}$C, which is a typical value in
winter. The effects of changing the undercooling temperature, $T_{t}$, are
qualitatively predictable and are not expected to change the occurrence of
different forms of the ice front morphology upon increasing the inclination
angle, $\beta$, as well as its physical explanation (by the boundary layer
theory based model). The only effects that $T_{t}$ will bring are that, 1) the
ice thickness at the equilibrium state of the system is thicker (thinner) with
decreasing (increasing) $T_{t}$; 2) the corresponding ice front morphology is
similar in the shape compared with that of the current study but will be
different in the extension and the local curvature.

We have conducted several typical sets of experiments (we only perform the
typical cases of the experiments and did not conduct every experiment
corresponding to the simulations because the experiments are time consuming):
1) for investigation of the effects of different heating condition, we have
performed $T_{b}=6^{\circ}$C, $10^{\circ}$C and $12^{\circ}$C, at the system
inclination angle of $\beta=90$; 2) for the investigation of the influence of
the cell inclination, we have done several different inclination angles,
$\beta$ (unit: deg), i.e., $\beta=0$, 50, 90, 180, at the hot plate
temperature $T_{b}=10^{\circ}$C;

In the experiments, the ice forms on the cold plate and grows in thickness
until the system reaches a statistical equilibrium state. The statistical
equilibrium state is reached when the standard deviation of the ice thickness
time series is less than 0.5%.

Changing the heating condition, $T_{b}$ can adjust the intensity of the two
counterrotating convective rolls. Recall that the water density reaches the
maximum, $\rho_{c}$, at $T_{c}$ (which is approximately $4^{\circ}$). Then
depending on the value of the hot plate temperature, $T_{b}$, the flow
configuration can be classified into two regimes:

1) When $4^{\circ}\text{C}<T_{b}\leq 8^{\circ}\text{C}$, the anticlockwise
convective roll (with the temperature ranging from $T_{c}$ to $T_{b}$ and the
corresponding thermal driving $T_{b}-T_{c}\leq 4$K) has weaker intensity than
that of the clockwise (with the temperature ranging from $T_{0}$ to $T_{c}$
and the corresponding thermal driving $T_{c}-T_{0}\approx 4$K). In this
regime, the whole ice front is shielded by the colder clockwise roll and is
away from the influence of the hot plumes from the hot plate. So the ice front
is flat with a tilting angle.

2) When $T_{b}>8^{\circ}\text{C}$, the anticlockwise convective roll (with the
temperature ranging from $T_{c}$ to $T_{b}$ and the corresponding thermal
driving $T_{b}-T_{c}>4$K) has stronger intensity than that of the clockwise
(with the temperature ranging from $T_{0}$ to $T_{c}$ and the corresponding
thermal driving $T_{c}-T_{0}\approx 4$K). So the anticlockwise roll penetrates
the clockwise roll and finally affects the ice front morphology, which
represents thinner in thickness with higher $T_{b}$ but similar shape.

With the knowledge of these, we later limited our investigation to the case
with $T_{b}=10^{\circ}$C and performed more experimental investigation under
different inclination angles (see Fig. S2).

Figure S2: More results of the experiments ((a), (c), (e)) compared with the
simulations (with properly considering the water density anomaly) ((b), (d),
(f)). Parameter chosen: (a) and (b): $T_{b}=10^{\circ}$C and $\beta=0$; (c)
and (d): $T_{b}=10^{\circ}$C and $\beta=50$; (e) and (f): $T_{b}=10^{\circ}$C
and $\beta=180$. The unit for $\beta$ is deg. The blue dashed line in (a),
(c), and (e) shows the ice front position which is guide to the eye. (b), (d),
and (f) are the temperature field overlapped with isotherms as well as the
velocity vectors. The ice front is represented by the black thick line and the
$T_{c}$-isotherm is represented by the red thick line. The blue shaded area is
the ice layer.

Fig. S2 reports more experimental results with the parameter
$T_{b}=10^{\circ}$C and $T_{t}=-10^{\circ}$C, with inclination angle,
$\beta=0$(panels (a) and (b)), $\beta=50$(panels (c) and (d)), and
$\beta=180$(panels (e) and (f)). Fig. S2(a), (c), and (e) are from
experiments, and Fig. S2(b), (d), and (f) are from the corresponding
simulations.

From Fig. S2, we can conclude that the results from the experiments and
numerical simulations under different inclination angles agree well with each
other.

##  Appendix B Numerical Methods

The simulations are performed by means of the CH4-PROJECT code Calzavarini
(2019), which adopts a Lattice-Boltzmann algorithm for the description of
fluid and temperature dynamics, and an enthalpy method for the ice evolution.
The method has been validated against experiments Wang  _et al._ (2021), while
the code has been intensively tested in Rabbanipour Esfahani  _et al._ (2018).

In the simulation, we use the Boussinesq approximation, which means that the
density is regarded as a constant value except for that in the buoyancy term
in the momentum equation. The non-monotonic relationship of the water density
as a function of the temperature is
$\rho_{\text{w}}=\rho_{c}(1-\alpha^{*}|T_{b}-4|^{q})$ Gebhart and Mollendorf
(1977), with $\rho_{c}=999.972~{}kg/m^{3}$ the maximum density corresponding
to $T_{c}\approx 4^{\circ}$C, $\alpha^{*}=9.30\times 10^{-6}(K^{-q})$, and
$q=1.895$. All the physical properties of water and ice phase, except for
$\rho_{w}$ in the buoyancy term are evaluated at the mean temperatures in each
phase which are $(T_{\text{b}}+0)/2$ and $(T_{\text{t}}+0)/2$, respectively.
In the simulations we neglect the microscopic physics leading to kinetic
undercooling, Gibbs-Thomson effect and the anisotropic growth/melting Dash
_et al._ (2006). Furthermore, we assume the ice and water density remain the
same, i.e., $\rho_{I}=\rho_{w}$, to satisfy the incompressibility of the flow
in the simulation. The relevant equations that govern the fluid flow and the
boundary conditions are reported here below (the energy equation will be
presented later),

$\begin{split}&\vec{\nabla}\cdot\vec{u}=0,\\\ &\frac{\partial\vec{u}}{\partial
t}+\vec{u}\cdot\vec{\nabla}\vec{u}=-\frac{\vec{\nabla}p}{\rho_{\text{c}}}+\nu_{w}\nabla^{2}\vec{u}+\alpha^{*}g|T-4|^{q}\hat{e}_{z},\end{split}$
(S1)

where $\vec{u}(x,y,t)$, $p(x,y,t)$ are fluid velocity, pressure, respectively
and we denote with $x$ as the horizontal direction and with $y$ the vertical
direction; $\nu_{w}$, $\rho$, $g$ are the kinematic viscosity of water, the
density, and the acceleration of gravity, respectively. $\hat{e}_{z}$ is the
unit vector pointing in the direction opposite to that of gravity.

Figure S3: Schematic description of the problem considered. The blue shaded
region corresponds to the solid phase (ice) and the red shaded region to the
liquid phase (water). $T_{0}$ is the melting temperature. The top plate is
cooled at constant temperature $T_{t}$, and the bottom plate is heated at the
constant temperature $T_{b}$. The ice-water interface is fixed at $T_{0}$ and
is of no slip condition.The two side walls are of adiabatic conditions.

The boundary conditions corresponding to the governing equations above are
isothermal at the top and bottom plates, no-slip at the bottom plate and at
the ice-water interface, adiabatic at the lateral boundaries, and no-slip and
freezing (namely, Stefan condition Bodenschatz  _et al._ (2000)) at the phase-
changing interface (see also figure S3). The boundary conditions read:

$\begin{split}&T(x,0,t)=T_{\text{b}},\\\ &T(x,1,t)=T_{\text{t}},\\\
&\vec{u}(x,0,t)=0,\\\ &\vec{u}(x,1-h_{\text{i}}(x,t),t)=0,\\\
&\frac{\partial{T(x,y,t)}}{\partial x}|_{x=0}=0,\\\
&\frac{\partial{T(x,y,t)}}{\partial x}|_{x=L_{x}/H}=0,\\\
&L\rho_{\text{i}}V_{n}=\vec{n}\cdot\vec{q}_{\text{w}}-\vec{n}\cdot\vec{q}_{\text{i}},\end{split}$
(S2)

where $L$ is the latent heat, $V_{n}$ is the normal speed of the ice-water
interface, and $h_{i}(x,t)$ the dimensionless ice layer thickness at the
position $x$, $\vec{q}$ is the heat flux vector, $\vec{n}$ is a unit normal at
the ice-water interface pointing into the liquid. The subscripts I and w refer
to the ice and the water, respectively. The heat flux reads
$q_{\text{i}}=-k_{\text{i}}\vec{\nabla}T_{\text{i}}$ and
$q_{\text{w}}=-k_{\text{w}}\vec{\nabla}T_{\text{w}}$.

The boundary condition at the ice-water interface requires particular care due
to its time and space dependent character. So an useful method is to separate
the total enthalpy $h$ into sensible heat and latent heat Faden  _et al._
(2019):

$\text{h}=\left\\{\begin{split}&L\phi_{w}+\text{C}_{\text{pi}}T,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}when~{}T<T_{0},\\\
&L\phi_{w}+\text{C}_{\text{pi}}T_{0},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}when~{}T=T_{0},\\\
&L\phi_{w}+\text{C}_{\text{pi}}T_{0}+\text{C}_{\text{pw}}(T-T_{0}),~{}~{}~{}~{}~{}when~{}T>T_{0}.\end{split}\right.$
(S3)

where $T_{0}$ is the phase change temperature ($T_{0}=0$), and
$\phi_{w}(x,y,t)$ is the liquid fraction in the system and the relation
between $h_{i}(x,t)$ and $\phi_{w}(x,y,t)$ is
$h_{i}(x,t)=1-\int_{0}^{1}\phi_{w}(x,y,t)\,dy$. In the ice phase,
$\phi_{w}=0$, and in the water phase, $\phi_{w}=1$, which leads to an
additional source term, $S_{1}$ (the formulation is shown later), from the
latent heat contribution at the ice-water interface in the energy conservation
equation.

We use the Lattice Boltzmann method (LBM) which is able to capture the
turbulent convective dynamics in the water phase and also describe the phase
change process at the ice-water interface. The basic principle and formulation
of the method are described in Succi (2001), while the specific applications
to phase-change have been presented in Huber  _et al._ (2008); Rabbanipour
Esfahani  _et al._ (2018). It is noteworthy that the key to accurately solve
such problems is to recover the diffusion term in the energy conservation
equation exactly and, similarly to Chen  _et al._ (2017), we implement the
correction when the investigated domain consists of heterogeneous media which
lead to another additional source term, $S_{2}$, in the energy conservation
equation. So the energy equation with consideration of two source terms
$S_{1}$ and $S_{2}$ reads

$\sigma(\rho\text{C}_{\text{p}})_{0}\frac{\partial T}{\partial
t}+\vec{\nabla}\cdot(\sigma(\rho\text{C}_{\text{p}})_{0}T\vec{u})=\vec{\nabla}\cdot(k\vec{\nabla}T)+S_{1}+S_{2},$
(S4)

where $T(x,y,t)$ is the temperature fields (all temperatures are measured in
Celsius), and $k$ is the thermal conductivity. When it is in water phase
$k=k_{w}$, $\text{C}_{\text{p}}=\text{C}_{\text{pw}}$, and when it is in ice
phase $k=k_{i}$, $\text{C}_{\text{p}}=\text{C}_{\text{pi}}$. The first source
term is $S_{1}=-L\rho\frac{\partial\phi_{w}}{\partial t}$ and the second
source term $S_{2}=-\sigma
k\vec{\nabla}T\vec{\nabla}\frac{1}{\sigma}-\frac{\rho\text{C}_{\text{p}}}{\sigma}T\vec{u}\vec{\nabla}\sigma$.
Here $\sigma=\frac{\rho\text{C}_{\text{p}}}{(\rho\text{C}_{\text{p}})_{0}}$ is
the ratio of heat capacitance (which is variable and depends on the type of
phase, i.e., ice or water) and $(\rho\text{C}_{\text{p}})_{0}$ is reference
heat capacitance which is taken as constant Chen  _et al._ (2017).

In this study we focus on the morphology of the ice front at the statistical
equilibrium state. In order to speed up the process of reaching an equilibrium
state, all the simulations are started with a thin layer of flat ice.

##  Appendix C Systematic investigation of various inclination angles

We have conducted a systematical exploration of how $\beta$ (unit:deg) affects
the ice morphology, ranging from $\beta=0$ to $\beta=180$ (here we use the
coordinate system attached to the cell). Fig. S4 reports the time sequences of
the temperature field under different inclination angles, $\beta$ with the
parameters $T_{t}=-10^{\circ}$C. We focus on the morphology of the ice front
at the statistical equilibrium state. In order to speed up the process of
reaching an equilibrium state, all the simulations are started with a thin
layer of flat ice.

The boundary conditions are no-slip for the velocity, adiabatic at the
sidewalls, and constant temperatures at the top and bottom plates. The initial
condition is still fluid at linear temperature profile from $T_{0}$ to
$T_{b}$. We assume thermophysical properties to be constant except for the
density in the buoyancy term. The real water density property near to the
density peak temperature, $T_{c}$, is well described with the equation
$\rho(T)=\rho_{c}(1-\alpha^{*}|T-T_{\text{c}}|^{q})$, with
$\rho_{c}=999.972~{}kg/m^{3}$ the maximum density corresponding to
$T_{c}\approx 4^{\circ}$C, $\alpha^{*}=9.30\times 10^{-6}(K^{-q})$, and
$q=1.895$ This form of equation gives the maximum density,
$\rho_{c}=999.972kg/m^{3}$, at the density peak temperature $T_{c}$ Gebhart
and Mollendorf (1977).

Figure S4: The systematic investigation of the influence of the inclination
angle on the ice front morphology. The snapshots at the internal, middle, and
outside circles are the temperature field 10min, 2.5hr, and 2.5d (when the
system has reached the statistical equilibrium state) after the start of the
simulations, respectively. The black thin lines are the isotherms. The red
dashed line shows the ice front, and the red thick line shows the
$T_{c}$-isotherm. The sketch in the middle shows the coordinate system
(attached to the system) and how the inclination angle, $\beta$ (unit: deg),
is defined.

Fig. S4 shows the temperature field (overlapped with isotherms) with varying
the system inclination angle, $\beta$ (unit: deg), from 0 (Rayleigh-Bénard
convection with heating from below and freshwater solidification from above)
to $180$ (flipped Rayleigh-Bénard convection with heating from above and
freshwater solidification from below), and here we use the coordinate system
attached to the cell. The snapshots at the internal, middle, and outside
circles are the temperature field 10min, 2.5hr, and 2.5d (when the system has
reached the statistical equilibrium state) after the start of the simulations,
respectively. In the RB (small $\beta$) and the flipped-RB (large $\beta$)
regimes, the ice front is heavily affected by the intensive convective flow,
and thus the ice front represents complex forms of morphology different from
those near the vertical convection cases.

##  Appendix D The boundary layer model and its input parameter

It has been mentioned that to locate the ice front profile, we need an input
from the simulation results, namely the ice thickness at the starting point of
the thermal boundary layer, $h_{i}(0)$. On top of this, the ice morphology can
be predicted.

Fig. S5 report the $h_{i}(0)$ as a function of the inclination angle, $\beta$.
It can been observed that $h_{i}(0)$ is roughly constant and around the level
of the averaged ice thickness $h_{i}$, which is shown in Fig. 3 of the main
paper.

Figure S5: $h_{i}(0)$ as a function of the inclination angle, $\beta$.
$h_{i}(0)$ is the ice thickness at $x=0$, which is the input parameter of the
boundary layer model.

##  Appendix E The buoyancy-intensity model

From the ice front morphology of the numerical simulation shown in
Figure.4(c)-(m) of the main paper, we can observe that along x-axis, the ice
thickness firstly increases and later decreases. We indicate with $\Delta
x_{0}$, the position where the ice thickness reaches the maximum (which is
here the a dimensionless position normalized by $H$).

The prescribed system inclination results in different levels of thermal
stratification, which induces modifications of the ice front morphology. When
$\beta=0$, there is a stably-stratified layer (from $T_{0}$ to $T_{c}$) on top
of the unstably-stratified layer (from $T_{c}$ to $T_{b}$). As the cell is
tilted with a small $\beta$, the convection is strong enough to squeeze the
stably-stratified layer to be closely attached to the ice front, so the ice
morphology is influenced by the one-roll convective flow. As $\beta$ increases
towards $90$, the inclination of the temperature gradient with respect to the
gravity is strong enough to break down the stratification and thus the stably-
stratified layer is set into motion in the form of a clockwise convective roll
(an upward cold water current) which competes with the initially anticlockwise
roll (downward warm water current). The shielding effect is prominent because
the ice thickness reaches a local maximum, and the hot plumes impacting on the
top part of the ice induces a local minimum of the ice thickness. The flow
motion of the original stably-stratified layer (from $T_{0}$ to $T_{c}$)
intensifies as $\beta$ increases but still when $\beta<90$. After $\beta=90$,
the stratification configuration flips over. The temperature difference of the
unstably-stratified layer is now $(T_{c}-T_{0})\approx 4$K. The unstably
stratified layer is directly in contact with the ice front.

From the above described flow structure dependence on the inclination, we can
deduce that there are two aspects of convection motion: 1) The intrinsic
aspect, which means that even without inclination the convection exists, and
this aspect of convective motion occurs in the initially unstably-stratified
layer, i.e., from $T_{c}$ to $T_{b}$ when $\beta<90$ and from $T_{0}$ to
$T_{c}$ when $\beta>90$. 2) The inclination-induced aspect, which means the
initially stably-stratified layer can be set into convective motion when the
inclination is introduced into the system. This aspect of convective motion
occurs in the initially stably-stratified layer, i.e., from $T_{0}$ to $T_{c}$
when $\beta<90$ and from $T_{c}$ to $T_{b}$ when $\beta>90$. That’s why in the
definition of the buoyancy intensity, $I_{1}(\beta)$ and $I_{2}(\beta)$, there
are two parts: the first part that includes $g_{x}$ results from the intrinsic
aspect, and the second part that includes $g_{y}$ results from the
inclination-induced aspect. In order to display $I_{1}(\beta)$ and
$I_{2}(\beta)$ in a more compact form, we introduce the Heaviside step
function, $\mathcal{H}$, to turn on or off the intrinsic aspect. So the final
definition of $I_{1}(\beta)$ and $I_{2}(\beta)$ is of this kind of form (see
Eqn (S5))

$\begin{split}&I_{1}(\beta)=g_{x}(1-\tfrac{\rho(T_{m})}{\rho_{c}})-g_{y}\mathcal{H}[\beta-90^{\circ}](1-\tfrac{\rho(T_{0})}{\rho_{c}});\\\
&I_{2}(\beta)=g_{x}(1-\tfrac{\rho(T_{m2})}{\rho_{c}})+g_{y}\mathcal{H}[90^{\circ}-\beta](1-\tfrac{\rho(T_{b})}{\rho_{c}}).\\\
\end{split}$ (S5)

The competition of the two convective rolls result in a local minimum heat
transfer rate where the ice thickness reaches the maximum, and we use the
buoyancy-intensity ratio to quantify this competition effect which give good
agreement with the trend we observed in the numerical simulations (see
Figure.4(b) of the main paper). This provides further evidence that the
competition of two convective rolls accounts for the form of ice front
morphology.

## References

  * Gebhart and Mollendorf (1977) B. Gebhart and J. C. Mollendorf, Deep Sea Research Part II Topical Studies in Oceanography 24, 831 (1977).
  * Calzavarini (2019) E. Calzavarini, Software Impacts 1, 100002 (2019).
  * Wang  _et al._ (2021) Z. Wang, E. Calzavarini, C. Sun, and F. Toschi, Proc. Natl. Acad. Sci. 118(10), e2012870118 (2021).
  * Rabbanipour Esfahani  _et al._ (2018) B. Rabbanipour Esfahani, S. C. Hirata, S. Berti, and E. Calzavarini, Phys. Rev. Fluids 3, 053501 (2018).
  * Dash  _et al._ (2006) J. Dash, A. Rempel, and J. Wettlaufer, Reviews of modern physics 78, 695 (2006).
  * Bodenschatz  _et al._ (2000) E. Bodenschatz, W. Pesch, and G. Ahlers, Annual review of fluid mechanics 32, 709 (2000).
  * Faden  _et al._ (2019) M. Faden, A. König-Haagen, and D. Brüggemann, Energies 12, 868 (2019).
  * Succi (2001) S. Succi,  _The Lattice-Boltzmann equation: for fluid dynamics and beyond_ (Oxford university press, 2001).
  * Huber  _et al._ (2008) C. Huber, A. Parmigiani, B. Chopard, M. Manga, and O. Bachmann, Int. J. Heat & Fluid Flow 29, 1469 (2008).
  * Chen  _et al._ (2017) S. Chen, Y. Y. Yan, and W. Gong, Int. J. Heat and Mass Transf. 107, 862 (2017).


# Polarization-controlled Brillouin scattering in elliptical optophononic
resonators

A. Rodriguez These authors contributed equally. Université Paris-Saclay, CNRS,
Centre de Nanosciences et de Nanotechnologies, 91120 Palaiseau, France  E.
Mehdi  Priya  Université Paris-Saclay, CNRS, Centre de Nanosciences et de
Nanotechnologies, 91120 Palaiseau, France  E. R. Cardozo de Oliveira
Université Paris-Saclay, CNRS, Centre de Nanosciences et de Nanotechnologies,
91120 Palaiseau, France  M. Esmann  Université Paris-Saclay, CNRS, Centre de
Nanosciences et de Nanotechnologies, 91120 Palaiseau, France  N. D.
Lanzillotti-Kimura<EMAIL_ADDRESS>Université Paris-Saclay,
CNRS, Centre de Nanosciences et de Nanotechnologies, 91120 Palaiseau, France

###### Abstract

The fast-growing development of optomechanical applications has motivated
advancements in Brillouin scattering research. In particular, the study of
high frequency acoustic phonons at the nanoscale is interesting due to large
range of interactions with other excitations in matter. However, standard
Brillouin spectroscopy schemes rely on fixed wavelength filtering, which
limits the usefulness for the study of tunable optophononic resonators. It has
been recently demonstrated that elliptical optophononic micropillar resonators
induce different energy-dependent polarization states for the Brillouin and
the elastic Rayleigh scattering, and that a polarization filtering setup could
be implemented to increase the contrast between the inelastic and elastic
scattering of the light. An optimal filtering configuration can be reached
when the polarization states of the laser and the Brillouin signal are
orthogonal from each other. In this work, we theoretically investigate the
parameters of such polarization-based filtering technique to enhance the
efficiency of Brillouin scattering detection. For the filtering optimization,
we explore the initial wavelength and polarization state of the incident
laser, as well as in the ellipticity of the micropillars, and reach an almost
optimal configuration for nearly background-free Brillouin detection. Our
findings are one step forward on the efficient detection of Brillouin
scattering in nanostructures for potential applications in fields such as
optomechanics and quantum communication.

##  1 Introduction

The engineering of nanomechanical systems has major impact in heat transport
management [1, 2], opto-acoustic cooling [3], quantum correlations mediated by
single phonons [4], and data storage [5, 6] among others. These advances
motivate the development of new techniques to investigate high frequency
phonons in the GHz range with high precision [7, 8, 9, 10, 11]. For instance,
Brillouin scattering, the inelastic scattering of light with acoustic phonons,
has experienced a renewed interest in the communities of nanophononics and
opto-nano-mechanics [12, 13].

The intrinsic properties of the materials govern the Brillouin scattering
selection rules and constrain the energy, direction, and polarization of the
Brillouin signal [14]. In addition, the overlap between optical and acoustic
fields determines the efficiency of the inelastic scattering process. In this
context, GaAs/AlAs optophononic systems allow to simultaneously control the
confinement of ultrahigh-frequency acoustic phonons and photons in the near-
infrared [15, 16, 17]. The overlap of the optical and acoustical fields in
such optophononic cavity leads to the enhancement of the Brillouin scattering
signals [18, 19]. Three-dimensional confinement of light and sound can be
achieved in semiconductor multilayered devices etched as micropillars [20,
21]. Such cavity systems receive attention as they can be, for example,
integrated with quantum dot to develop single photon sources [22, 23, 24]. To
detect the spontaneous Brillouin scattering signals, a particularly stringent
filtering of the elastically scattered laser is necessary. Therefore, new
techniques of tunable Brillouin spectroscopy are developed to allow the
efficient filtering of the scattered signal [25, 19, 20].

Here, we study elliptical micropillars, that enable the control of the
Brillouin scattering polarization selection rules [26]. Indeed, the anisotropy
of the micropillar cross-section lifts the degeneracy of the optical cavity
modes, in energy and polarization [27, 28, 29]. A wavelength-dependent
polarization rotation is then induced by the birefringence of the elliptical
microcavity [30, 31]. Contrary to standard polarization selection rules in
GaAs grown along the [001] crystalline direction, the Brillouin scattering
polarization in elliptical micropillars can be rotated differently than the
polarization of the excitation laser reflected by the optical cavity.

In a previous work [26], we reported the theoretical and experimental
demonstration of the manipulation of the Brillouin scattering polarization
with an excitation laser polarized along the linear diagonal polarization. In
this work, we theoretically introduce different parameters to control the
polarization of the Brillouin scattering.We demonstrate the possibility to
efficiently filter out the reflected laser by reaching an almost cross-
polarization configuration between the laser and the Brillouin signal.

##  2 Elliptical optophononic micropillar cavities

We theoretically study here an optophononic microcavity based on distributed
Bragg reflectors (DBRs), grown on a (001)-oriented GaAs substrate by
molecular-beam epitaxy. We use experimental data from the same sample as in
[26] to base our study on experimentally realistic conditions. A resonant
spacer of optical path length of $\lambda$/2 at a resonance wavelength of
$\lambda\sim 900\,$nm is enclosed by the two DBRs. The top (bottom) optical
DBR is formed by 25 (29) periods of
$\mathrm{Ga}_{\mathrm{0.9}}\mathrm{Al}_{\mathrm{0.1}}\mathrm{As}$/$\mathrm{Ga}_{\mathrm{0.05}}\mathrm{Al}_{\mathrm{0.95}}\mathrm{As}$
bilayers. Elliptical micropillars as schematised in Fig. 1a are fabricated
from this planar optophononic cavity by optical lithography followed by
Inductively Coupled Plasma Etching (ICP). A schematic of the elliptical cross-
section is shown in Fig. 1b where the major and minor axes are indicated (m
and n, respectively). The optical cavity presents two eigenmodes with
orthogonal linear polarization and two nondegenerate eigenenergies,
$\omega_{\mathrm{cav}}^{H}$ and $\omega_{\mathrm{cav}}^{V}$ (considering
$\hbar=1$), due to the anisotropy of the micropillar cross-section. We define
as well, for each mode, its corresponding wavelength
$\lambda_{\mathrm{cav}}^{H/V}$. These two optical modes correspond to the
horizontal ($H$) and vertical ($V$) polarizations.

When considering a linearly polarized incident laser beam,
$\ket{\psi_{\mathrm{in}}}$, the output polarization
$\ket{\psi_{\mathrm{out}}}$ of the laser reflected by the cavity is rotated as
a function of the laser wavelength $\lambda_{\mathrm{laser}}$ due to the
cavity birefringence [30, 32]. Both polarization states can be defined as
follows:

$\begin{split}&\ket{\psi_{\mathrm{in}}}=\dfrac{1}{\sqrt{|b_{\mathrm{in}}^{H}|^{2}+|b_{\mathrm{in}}^{V}|^{2}}}(b_{\mathrm{in}}^{H}\ket{H}+b_{\mathrm{in}}^{V}\ket{V}),\\\
&\ket{\psi_{\mathrm{out}}}=\dfrac{1}{\sqrt{|b_{\mathrm{out}}^{H}|^{2}+|b_{\mathrm{out}}^{V}|^{2}}}(b_{\mathrm{out}}^{H}\ket{H}+b_{\mathrm{out}}^{V}\ket{V}),\end{split}$
(1)

where $b_{\mathrm{in}}^{H/V}$ and $b_{\mathrm{out}}^{H/V}$ are the average
values of the polarization dependent input and output external field operators
[26]. These operators are defined with the input-output relation [33]:

$b_{\mathrm{out}}^{H/V}=b_{\mathrm{in}}^{H/V}+\sqrt{\kappa_{\mathrm{top},H/V}}a_{H/V}$
(2)

with $\kappa_{\mathrm{top},H/V}$ being the top DBR damping rate and $a_{H/V}$
the average value of the intracavity field operator. The reflected light,
represented by $b_{\mathrm{out}}^{H/V}$, is the result of the interference
between the light emerging from the cavity and the light directly reflected by
the top DBR. For a given input polarization $\ket{\psi_{\mathrm{in}}}$, the
response in polarization of the cavity can be calculated with the
polarization-dependent complex reflection coefficients $r_{H}$ and $r_{V}$
[34]:

$\begin{split}r_{H/V}=&\dfrac{b_{\mathrm{out}}^{H/V}}{b_{\mathrm{in}}^{H/V}}\\\
=&1-2\eta_{\mathrm{top},H/V}\dfrac{1}{1-\dfrac{2i(\omega_{\mathrm{in}}-\omega_{\mathrm{cav}}^{H/V})}{\kappa_{H/V}}},\end{split}$
(3)

where $\omega_{\mathrm{in}}$ and $\omega_{\mathrm{cav},H/V}$ are the energies
of the incident laser and the cavity modes, respectively. The output coupling
efficiency through the top DBR $\eta_{\mathrm{top},H/V}$ is defined as
$\eta_{\mathrm{top},H/V}=\kappa_{\mathrm{top},H/V}/\kappa_{H/V}$ with
$\kappa_{H/V}$ the total cavity damping rate including leakage through the top
and bottom DBRs, as well as the sidewall losses.

To theoretically study the polarization rotation induced by the elliptical
cavities, we use the experimental optical cavity parameters of a given
micropillar. A micropillar of ellipticity $e=\sqrt{m/n}-1=0.29$ is chosen from
the sample used in reference [26]. The two optical cavity modes can be then
characterized with reflectivity measurements as a function of the incident
laser wavelength as shown in Fig. 1c. To measure the reflectivity spectra of
the micropillar cavity, the polarization of the excitation laser is set to
match one of the cavity eigenmodes i.e., $\ket{\psi_{\mathrm{in}}}=\ket{H/V}$.
The experimental reflectivities (dots) of the $H$ mode (blue) and the $V$ mode
(red) are fitted (line) with the optical cavity parameters
($\lambda_{\mathrm{cav},H/V}$, $\eta_{\mathrm{top},H/V}$ and $\kappa_{H/V}$)
as the cavity reflectivity is directly deduced from $|r_{H/V}|^{2}$ for each
mode.

Figure 1: (a) Schematic of the micropillar vertical layer structure with two
distributed Bragg reflectors (blue) surrounding a resonant spacer (light
yellow). (b) Schematic of the micropillar cross-section with $m$ and $n$ the
major and minor axes of the elliptical structure. The two linear polarizations
$H$ and $V$ are the two fundamental optical eigenmodes of the micropillar. (c)
Measured (dots) and fitted (solid line) optical reflectivity of an elliptical
pillar with an ellipticity of e = 0.29. The blue (red) line corresponds to
$\ket{H}$ ($\ket{V}$).

##  3 Polarization rotation induced by the elliptical cavities

For an arbitrary linear input polarization state $\ket{\psi_{\mathrm{in}}}$,
we calculate the polarization rotation induced by the considered micropillar
cavity with Eq. (3) and plot it in the Poincaré sphere as a function of the
input laser wavelength (see colorscale) as shown in Fig. 2a. For each laser
wavelength, the Stokes parameters of the output state
$\ket{\psi_{\mathrm{out}}}$ are calculated and give the coordinates of the
polarization state in the Poincaré sphere. Far from the cavity resonance,
there is no polarization rotation, i.e.,
$\ket{\psi_{\mathrm{out}}}=\ket{\psi_{\mathrm{in}}}$, while
$\ket{\psi_{\mathrm{out}}}$ forms a closed path on the surface of the Poincaré
sphere when in resonance with the cavity, assuming a perfect mode matching
between the input beam and optical cavity. The trajectory of the output
polarization state as a function of the input laser wavelength depends
directly on the cavity parameters and the laser energy $\omega_{\mathrm{in}}$.

The simulations of the output polarization states are based on the Jones
matrices formalism [35, 26]. The cavity parameters deduced from the fit of the
cavity reflectivity spectra are used as input for the Jones matrices to
calculate the polarization rotation induced by the cavity.

The polarization of the Brillouin scattering field, inside the GaAs cavity
spacer, is calculated for every value of the excitation laser wavelength. Its
polarization amplitude is the same as the intracavity polarization state of
the input laser but at a different energy $\omega_{\mathrm{B}}$. The Brillouin
scattering polarization state at the output of the cavity,
$\ket{\psi_{\mathrm{B}}}_{\mathrm{B/AS}}$, is thus defined as:

$\ket{\psi_{\mathrm{B}}}_{\mathrm{B/AS}}=\dfrac{1}{\sqrt{|b_{\mathrm{B}}^{H}|^{2}+|b_{\mathrm{B}}^{V}|^{2}}}(b_{\mathrm{B}}^{H}\ket{H}+b_{\mathrm{B}}^{V}\ket{V}),$
(4)

with $b_{\mathrm{B}}^{H}$ the average value of the Brillouin scattering
operator defined as follow:

$\begin{split}&b_{\mathrm{B}}^{H}\propto\\\
&\dfrac{1}{1-\dfrac{2i(\omega_{\mathrm{B}}-\omega_{\mathrm{cav}}^{H/V})}{\kappa_{H/V}}}\cdot\dfrac{1}{1-\dfrac{2i(\omega_{\mathrm{in}}-\omega_{\mathrm{cav}}^{H/V})}{\kappa_{H/V}}}\cdot
b_{\mathrm{in}}^{H/V}\end{split}$ (5)

The first term depends on the frequency of the Brillouin signal
$\omega_{\mathrm{B}}$. In our case, we consider only the first harmonic at
$\pm 18\,$GHz which corresponds to the confined acoustic mode in the spacer.
The second term depends on the incoming laser frequency $\omega_{\mathrm{in}}$
and its polarization with $b_{\mathrm{in}}^{H/V}$. The Brillouin scattering
polarization state at the output of the cavity is then calculated as a
function of the incoming laser wavelength for the Stokes and the anti-Stokes
signals at $18\,$GHz and plotted in the Poincaré sphere in the Fig. 2b. The
behaviour of the Brillouin scattering polarization is different from the one
of the light reflected by the cavity. The Brillouin signal can therefore be
separated from the reflected laser by polarization filtering.

From the Stokes parameters describing the reflected laser and the Brillouin
scattering polarization states, the angle $\theta$ between the polarization
states $\ket{\psi_{\mathrm{out}}}$ and
$\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$ in the Poincaré sphere is calculated
for the Stokes and anti-Stokes signals as follows:

$\theta=2\mathrm{arcsin}(\dfrac{\Delta}{2}),$ (6)

with
$\Delta=\sqrt{|s_{HV}^{\mathrm{B}}-s_{HV}^{\mathrm{out}}|^{2}+|s_{DA}^{\mathrm{B}}-s_{DA}^{\mathrm{out}}|^{2}+|s_{RL}^{\mathrm{B}}-s_{RL}^{\mathrm{out}}|^{2}}$.
The angle $\theta$ is plotted in Fig. 2c and 2d as a function of the
wavelength of the input laser for the Stokes and anti-Stokes signals
respectively. The ideal filtering condition of the reflected laser corresponds
to the two polarization states $\ket{\psi_{\mathrm{out}}}$ and
$\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$ orthogonal to each other, i.e.,
giving two opposite points in the Poincaré sphere. Under these conditions,
with an angle of $\theta=180^{\circ}$ between the two polarization states in
the Poincaré sphere, the reflected laser can be entirely filtered while
getting all of the Brillouin signal by selecting only the polarization along
$\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$. A maximal filtering of the
reflected laser can however be achieved when the Brillouin signal polarization
state is the furthest from the reflected laser in the Poincaré sphere,
corresponding to the maximal angle $\theta_{\mathrm{max}}$. In Fig. 2c and 2d,
the maximal angles reached for the Stokes and anti-Stokes signals,
$\theta_{\mathrm{max}}^{\mathrm{S}}=85^{\circ}$ and
$\theta_{\mathrm{max}}^{\mathrm{AS}}=102^{\circ}$, are indicated with the
corresponding wavelength $\lambda_{\mathrm{opt}}^{\mathrm{S}}$ and
$\lambda_{\mathrm{opt}}^{\mathrm{AS}}$. In this situation, the polarization of
the Stokes and anti-Stokes signal at $\pm 18\,$GHz are not orthogonal to the
reflected laser polarization as the maximal angle does not reach the value of
$180^{\circ}$. The polarization states $\ket{\psi_{\mathrm{out}}}$ and
$\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$ at the wavelength
$\lambda_{\mathrm{opt}}^{\mathrm{S}}$ and
$\lambda_{\mathrm{opt}}^{\mathrm{AS}}$ are now compared in the Poincaré sphere
in Fig. 2e and 2f, where the input polarization state
$\ket{\psi_{\mathrm{in}}}$ is plotted as well. It allows the direct comparison
of the polarization states $\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$ and
$\ket{\psi_{\mathrm{out}}}$ as they are not opposite to each other in the
sphere. A cross-polarization scheme, adjusted to the reflected laser
polarization, would allow to completely eliminated the laser. However, the
Brillouin signal will be highly attenuated in this situation, as its
polarization is not orthogonal to the laser.

For this given input polarization state $\ket{\psi_{\mathrm{in}}}$, the
optimal condition reached between the polarization states
$\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$ and $\ket{\psi_{\mathrm{out}}}$ is
far from ideal as the two polarization states are not opposite to each other
in the Poincaré sphere. However, these two polarization states depend by
definition (see Eq.(2) and Eq.(5)) on the average, i.e., expectation values,
of the polarization dependent external input field operators
$b_{\mathrm{in}}^{H/V}$, i.e., they depend on the input polarization state
$\ket{\psi_{\mathrm{in}}}$. Thus the angle $\theta$ can be furthermore
optimised by changing the input polarization state.

Figure 2: (a) Simulated polarization rotation of the polarization state
$\ket{\psi_{\mathrm{out}}}$ induced by the cavity birefringence displayed in
the Poincaré sphere as a function of the input laser wavelength indicated by
the colorscale. (b) Simulated polarization states of the Stokes and anti-
Stokes signal of the Brillouin scattering at $18\,$GHz,
$\ket{\psi_{\mathrm{B}}}_{\mathrm{S}}$ and
$\ket{\psi_{\mathrm{B}}}_{\mathrm{AS}}$, as a function of the input laser
wavelength. (c) and (d) Angle in the Poincaré sphere between the Stokes (resp.
anti-Stokes) signal polarization state and the laser polarization rotated by
the cavity as a function of the laser wavelength. The optimal wavelength
$\lambda_{\mathrm{opt}}^{\mathrm{S}}$ (resp.
$\lambda_{\mathrm{opt}}^{\mathrm{AS}}$) corresponds to the maximal angle
between $\ket{\psi_{\mathrm{B}}}_{\mathrm{S}}$ (resp.
$\ket{\psi_{\mathrm{B}}}_{\mathrm{AS}}$) and $\ket{\psi_{\mathrm{out}}}$. (e)
and (f) Representation in the Poincaré sphere of the polarization states
$\ket{\psi_{\mathrm{in}}}$, $\ket{\psi_{\mathrm{B}}}_{\mathrm{S}}$ (resp.
$\ket{\psi_{\mathrm{B}}}_{\mathrm{AS}}$) and $\ket{\psi_{\mathrm{out}}}$ at
the optimal wavelength $\lambda_{\mathrm{opt}}^{\mathrm{S}}$ (resp.
$\lambda_{\mathrm{opt}}^{\mathrm{AS}}$).

##  4 Optimization of the filtering

Now, we study the influence of the input polarization state
$\ket{\psi_{\mathrm{in}}}$ on the filtering of the laser while the micropillar
characteristics ($\lambda_{\mathrm{cav},H/V}$, $\eta_{\mathrm{top},H/V}$ and
$\kappa_{H/V}$) remain the same. The whole Poincaré sphere is covered by
varying $\ket{\psi_{\mathrm{in}}}$. For each given input polarization state,
the maximal angle $\theta_{\mathrm{max}}$ reached between the two polarization
states $\ket{\psi_{\mathrm{out}}}$ and
$\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$ in the Poincaré sphere is
calculated. In Fig. 3, this maximal angle is plotted on the Poincaré sphere
for the Stokes (Fig. 3.a) and the anti-Stokes (Fig. 3.b) signals. The
coordinates of each point on the Poincaré sphere correspond to the Stokes
parameters of the input polarization state $\ket{\psi_{\mathrm{in}}}$. The
value of the maximal angle $\theta_{\mathrm{max}}$ is then given by the
colorscale. It is important to note that $\theta_{\mathrm{max}}$ depends both
on the incident polarization $\ket{\psi_{\mathrm{in}}}$ and on the laser
wavelength $\lambda_{\mathrm{laser}}$. Therefore, in Fig. 3, the plotted value
of $\theta_{\mathrm{max}}$ is calculated for the laser wavelength which gives
the best filtering, $\lambda_{\mathrm{opt}}^{\mathrm{S}}$ and
$\lambda_{\mathrm{opt}}^{\mathrm{AS}}$, respectively. Hence, for each point,
i.e., each $\ket{\psi_{\mathrm{in}}}$, $\theta_{\mathrm{max}}$ is reached for
a different value of the input laser wavelength.

The minimal angle $\theta=0$ implies that
$\ket{\psi_{\mathrm{out}}}=\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$. It is is
reached when the input polarization state is aligned with the cavity
eigenmodes $\ket{H}$ or $\ket{V}$, as there is no polarization rotation
induced by the cavity in these cases. Therefore, if
$\ket{\psi_{\mathrm{in}}}=\ket{H/V}$,
$\ket{\psi_{\mathrm{in}}}=\ket{\psi_{\mathrm{out}}}=\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$.

The maximal angle $\theta_{\mathrm{max}}$ reached between the reflected laser
polarization state and the Brillouin signals is $105^{\circ}$. It is obtained
for incident polarizations $\ket{\psi_{\mathrm{in}}}$ with coordinates in the
Poincaré sphere corresponding to an annular shape around the HV axis and
parallel to the DRAL meridian (in pink in Fig. 3). This maximal angle does not
correspond to two opposite points in the Poincaré sphere. Thus, a cross-
polarization scheme, with the reflected laser polarization as reference, would
completely filter out the laser but attenuate as well the Brillouin signal.
For this given set of parameters ($\lambda_{\mathrm{cav},H/V}$,
$\eta_{\mathrm{top},H/V}$ and $\kappa_{H/V}$), there is no input polarization
state $\ket{\psi_{\mathrm{in}}}$ for which the states
$\ket{\psi_{\mathrm{out}}}$ and $\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$ are
opposite to each other in the Poincaré sphere.

Figure 3: Maximal angle $\theta_{\mathrm{max}}$ in the Poincaré sphere between
the polarization states of the reflected laser $\ket{\psi_{\mathrm{out}}}$ and
the Brillouin Stokes (a) and anti-Stokes (b) signals,
$\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$. The value of the angle is indicated
by the colorscale and is plotted as a function of the input polarization state
$\ket{\psi_{\mathrm{in}}}$.

The parameters studied so far ($\lambda_{\mathrm{laser}}$ and
$\ket{\psi_{\mathrm{in}}}$) can experimentally be tuned for a given
micropillar to reach the maximal angle between the two polarization states
$\ket{\psi_{\mathrm{out}}}$ and $\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$.
However, the orthogonality between those two polarization states cannot be
reached in this case. We can overcome this limitation by varying the
birefringence of the micropillar cross-section. Indeed, the splitting
$\Delta\lambda$ between the two cavity modes, which has a determinant role in
the rotation of polarization, is directly controlled by the ellipticity of the
micropillar [36]. Experimentally, the birefringence is monitored by etching
arrays of micropillars of varying ellipticity [26]. In the following, we study
the effect of the ellipticity on the output polarization states by modifying
$\Delta\lambda$ in the numerical simulations.

The maximal angle $\theta_{\mathrm{max}}$ in the Poincaré sphere, between the
two polarization states $\ket{\psi_{\mathrm{out}}}$ and
$\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$, depends on how much the two optical
cavity modes overlap in energy, which directly depends on $\Delta\lambda$, the
separation between the central energies of the two optical cavity modes. In
addition, the overlap between the two modes also depends on the full width at
half maximum defined by the total cavity damping rate $\kappa_{H/V}$. A large
overlap implies a small cavity-induced polarization rotation, and thus a small
$\theta$. In contrast, a small mode overlap results in a more important
rotation of polarization. In Fig. 4.a and b, we plot the maximal angle
$\theta_{\mathrm{max}}$ reached, associated to a given
$\ket{\psi_{\mathrm{in}}}$ and $\lambda_{\mathrm{laser}}$, as a function of
the cavity splitting between the two cavity modes $\Delta\lambda$ and for
different values of the total cavity damping rate $\kappa$, where
$\kappa_{H}=\kappa_{V}=\kappa$. The case where $\kappa=110\,\mu$eV (72 pm) and
$\Delta\lambda=0.06\,$nm corresponds to the situation studied above and is
indicated by a black square in Fig. 4. For each set of $\Delta\lambda$ and
$\kappa$, the maximal value $\theta_{\mathrm{max}}$ is calculated. We observe
that $\theta_{\mathrm{max}}$ increases with $\Delta\lambda$ and inversely with
$\kappa$. It reaches a maximum of $\theta_{\mathrm{max}}=170^{\circ}$ for
$\kappa=110\,\mu$eV (72 pm). For low values of cavity damping rate $\kappa$,
the spectral separation of the two modes increases rapidly with
$\Delta\lambda$. Consequently, we will reach $\theta_{\mathrm{max}}$ and
optimal conditions for a smaller ellipticity than in the case of a cavity with
high damping rate.

Figure 4:  (a) and (b) Maximal angle $\theta_{\mathrm{max}}$ between
$\ket{\psi_{\mathrm{out}}}$ and $\ket{\psi_{\mathrm{B}}}_{\mathrm{S/AS}}$
reached in the Poincaré sphere for a given $\ket{\psi_{\mathrm{in}}}$ for
different values of $\kappa$, the total cavity damping rate.

Experimentally, an optimal configuration can be achieved for a given
micropillar by tuning $\omega_{\mathrm{in}}$ and $\ket{\Psi_{\mathrm{in}}}$.
In addition, the ellipticity of the micropillar cross-section can be designed
through the fabrication process to control $\Delta\lambda$ and thus maximize
the angle $\theta_{\mathrm{max}}$. The total cavity damping rate
$\kappa_{H/V}$, on which the optimization of the filtering process depends as
well, can also be designed with the fabrication process. Indeed, the number of
layers in the DBRs controls the value of $\kappa_{H/V}$, yet limited by the
micropillar fabrication. Moreover, the spatial mode matching between the
incoming laser beam and a micropillar can be experimentally limited by the
ellipticity of the micropillar, which degrades the coupling between the input
light and the cavity. In such a situation, the efficiency and the intensity of
the Brillouin scattering is compromised. Thus, finding the balance between an
efficient mode matching between the optical cavity and the input laser, and
the ellipticity of the cross-section is necessary to implement this
polarization-based Brillouin spectroscopy technique.

##  5 Conclusion

We have theoretically studied the polarization control of elliptical
optophononic micropillar resonators to optimize the Brillouin scattering
detection. We demonstrated that we can maximize the angle between the
Brillouin signal and the excitation laser polarization states by independently
manipulating them. The reflected laser filtering can therefore be implemented
with a cross-polarization scheme, while attenuating the Brillouin scattered
signal. However, the numerical simulations presented in this work are
restricted as they only take into account the rotation of polarization and not
the intensity of the Brillouin scattering. It is necessary to consider the
enhancement of the Brillouin signal by the optical cavity. This study can be
experimentally applied to micropillars with different ellipticities, and be
performed using a similar setup as presented in [26], adding the control of
the incident polarization.

This Brillouin spectroscopy technique can be extended to other polarization-
dependent structures, such as nanofibers [37]. Therefore, our findings are
building blocks for efficient Brillouin detection in tunable optophononic
resonators for future opto-mechanical applications.

## Acknowledgements

The authors acknowledge funding from European Research Council Consolidator
Grant No.101045089 (T-Recs). This work was supported by the European
Commission in the form of the H2020 FET Proactive project No. 824140 (TOCHA),
and the French RENATECH network.

## References

  * [1] O. Florez, G. Arregui, M. Albrechtsen, R. C. Ng, J. Gomis-Bresco, S. Stobbe, C. M. Sotomayor-Torres, and P. D. García.  Engineering nanoscale hypersonic phonon transport.  Nature Nanotechnology, 17:947–951, 2022.
  * [2] Oliver Braun, Roman Furrer, Pascal Butti, Kishan Thodkar, Ivan Shorubalko, Ilaria Zardo, Michel Calame, and Mickael L. Perrin.  Spatially mapping thermal transport in graphene by an opto-thermal method.  npj 2D Materials and Applications, 6(1):6, 2022.
  * [3] Laura Blázquez Martínez, Philipp Wiedemann, Changlong Zhu, Andreas Geilen, and Birgit Stiller.  Optoacoustic cooling of traveling hypersound waves.  Physical Review Letters, 132:023603, 2024.
  * [4] Mitchell D. Anderson, Santiago Tarrago Velez, Kilian Seibold, Hugo Flayac, Vincenzo Savona, Nicolas Sangouard, and Christophe Galland.  Two-color pump-probe measurement of photonic quantum correlations mediated by a single phonon.  Physical Review Letters, 120:233601, 2018.
  * [5] Birgit Stiller, Moritz Merklein, Christian Wolff, Khu Vu, Pan Ma, Stephen J. Madden, and Benjamin J. Eggleton.  Coherently refreshing hypersonic phonons for light storage.  Optica, 7(5):492, 2020.
  * [6] Moritz Merklein, Birgit Stiller, Khu Vu, Pan Ma, Stephen J. Madden, and Benjamin J. Eggleton.  On-chip broadband nonreciprocal light storage.  Nanophotonics, 10:75–82, 2020.
  * [7] Sebastian Volz, Jose Ordonez-Miranda, Andrey Shchepetov, Mika Prunnila, Jouni Ahopelto, Thomas Pezeril, Gwenaelle Vaudel, Vitaly Gusev, Pascal Ruello, Eva M. Weig, Martin Schubert, Mike Hettich, Martin Grossman, Thomas Dekorsy, Francesc Alzina, Bartlomiej Graczykowski, Emigdio Chavez-Angel, J. Sebastian Reparaz, Markus R. Wagner, Clivia M. Sotomayor-Torres, Shiyun Xiong, Sanghamitra Neogi, and Davide Donadio.  Nanophononics: State of the art and perspectives.  The European Physical Journal B, 89(1):15, 2016.
  * [8] Priya, E. R. Cardozo de Oliveira, and N. D. Lanzillotti-Kimura.  Perspectives on high-frequency nanomechanics, nanoacoustics, and nanophononics.  Appl. Phys. Lett., 122:140501, 2023.
  * [9] E. R. Cardozo de Oliveira, C. Xiang, M. Esmann, N. Lopez Abdala, M. C. Fuertes, A. Bruchhausen, H. Pastoriza, B. Perrin, G. J. A. A. Soler-Illia, and N. D. Lanzillotti-Kimura.  Probing gigahertz coherent acoustic phonons in Tio2 mesoporous thin films.  Photoacoustics, 30:100472, 2023.
  * [10] O. Ortiz, F. Pastier, A. Rodriguez, Priya, A. Lemaitre, C. Gomez-Carbonell, I. Sagnes, A. Harouri, P. Senellart, V. Giesz, M. Esmann, and N. D. Lanzillotti-Kimura.  Fiber-integrated microcavities for efficient generation of coherent acoustic phonons.  Appl. Phys. Lett., 117:183102, 2020.
  * [11] Diego H.O. Machado, Antonio Crespo-Poveda, Alexander S. Kuznetsov, Klaus Biermann, Luis V.A. Scalvi, and Paulo V. Santos.  Generation and Propagation of Superhigh-Frequency Bulk Acoustic Waves in Ga As.  Physical Review Applied, 12(4):044013, 2019.
  * [12] Benjamin J. Eggleton, Christopher G. Poulton, Peter T. Rakich, Michael. J. Steel, and Gaurav Bahl.  Brillouin integrated photonics.  Nature Photonics, 13:664–677, 2019.
  * [13] Fariborz Kargar and Alexander A. Balandin.  Advances in Brillouin–Mandelstam light-scattering spectroscopy.  Nature Photonics, 15(10):720–731, 2021.
  * [14] Peter Y. Yu and Manuel Cardona.  Fundamentals of Semiconductors.  Springer Berlin Heidelberg, 2010.
  * [15] O. Ortiz, Priya, A. Rodriguez, A. Lemaître, M. Esmann, and N. D. Lanzillotti-Kimura.  Topological optical and phononic interface mode by simultaneous band inversion.  Optica, 8(5):598–605, 2021.
  * [16] A. Fainstein, N. D. Lanzillotti-Kimura, B. Jusserand, and B. Perrin.  Strong optical-mechanical coupling in a vertical GaAs/AlAs microcavity for subterahertz phonons and near-infrared light.  Phys. Rev. Lett., 110:037403, 2013.
  * [17] S. Anguiano, A. E. Bruchhausen, B. Jusserand, I. Favero, F. R. Lamberti, L. Lanco, I. Sagnes, A. Lemaître, N. D. Lanzillotti-Kimura, P. Senellart, and A. Fainstein.  Micropillar Resonators for Optomechanics in the Extremely High 19–95-GHz Frequency Range.  Physical Review Letters, 118(26):263901, 2017.
  * [18] N. D. Lanzillotti-Kimura, A. Fainstein, B. Jusserand, and A. Lemaître.  Resonant Raman scattering of nanocavity-confined acoustic phonons.  Physical Review B, 79:035404, 2009.
  * [19] A. Rodriguez, P. Priya, O. Ortiz, P. Senellart, C. Gomez-Carbonell, A. Lemaître, M. Esmann, and N. D. Lanzillotti-Kimura.  Fiber-based angular filtering for high-resolution Brillouin spectroscopy in the 20-300 GHz frequency range.  Optics Express, 29:2637, 2021.
  * [20] M. Esmann, F. R. Lamberti, A. Harouri, L. Lanco, I. Sagnes, I. Favero, G. Aubin, C. Gomez-Carbonell, A. Lemaître, O. Krebs, P. Senellart, and N. D. Lanzillotti-Kimura.  Brillouin scattering in hybrid optophononic Bragg micropillar resonators at 300 GHz.  Optica, 6:854, 2019.
  * [21] F. R. Lamberti, Q. Yao, L. Lanco, D. T. Nguyen, M. Esmann, A. Fainstein, P. Sesin, S. Anguiano, V. Villafañe, A. Bruchhausen, P. Senellart, I. Favero, and N. D. Lanzillotti-Kimura.  Optomechanical properties of GaAs/AlAs micropillar resonators operating in the 18 GHz range.  Optics Express, 25:24437–24447, 2017.
  * [22] N. Somaschi, V. Giesz, L. De Santis, J. C. Loredo, M. P. Almeida, G. Hornecker, S. L. Portalupi, T. Grange, C. Antón, J. Demory, C. Gómez, I. Sagnes, N. D. Lanzillotti-Kimura, A. Lemaítre, A. Auffeves, A. G. White, L. Lanco, and P. Senellart.  Near-optimal single photon sources in the solid state.  Nature Photon., 10:340–345, 2016.
  * [23] Hui Wang, Yu-Ming He, T.-H. Chung, Hai Hu, Ying Yu, Si Chen, Xing Ding, M.-C. Chen, Jian Qin, Xiaoxia Yang, Run-Ze Liu, Z.-C. Duan, J.-P. Li, S. Gerhardt, K. Winkler, J. Jurkat, Lin-Jun Wang, Niels Gregersen, Yong-Heng Huo, Qing Dai, Siyuan Yu, Sven Höfling, Chao-Yang Lu, and Jian-Wei Pan.  Towards optimal single-photon sources from polarized microcavities.  Nature Photonics, 13(11):770–775, 2019.
  * [24] Morten P. Bakker, Ajit V. Barve, Thomas Ruytenberg, Wolfgang Löffler, Larry A. Coldren, Dirk Bouwmeester, and Martin P. van Exter.  Polarization degenerate solid-state cavity quantum electrodynamics.  Phys. Rev. B, 91:115319, 2015.
  * [25] Guillermo Rozas, Bernard Jusserand, and Alejandro Fainstein.  Fabry-Pérot-multichannel spectrometer tandem for ultra-high resolution Raman spectroscopy.  Rev. Sci. Instrum., 85:013103, 2014.
  * [26] A. Rodriguez, P. Priya, E. R. Cardozo De Oliveira, A. Harouri, I. Sagnes, F. Pastier, L. Le Gratiet, M. Morassi, A. Lemaître, L. Lanco, M. Esmann, and N. D. Lanzillotti-Kimura.  Brillouin scattering selection rules in polarization-sensitive photonic resonators.  ACS Photonics, 10:1687–1693, 2023.
  * [27] M. T. Rakher, N. G. Stoltz, L. A. Coldren, P. M. Petroff, and D. Bouwmeester.  Externally mode-matched cavity quantum electrodynamics with charge-tunable quantum dots.  Phys. Rev. Lett., 102:097403, 2009.
  * [28] B. Gayral, J. M. Gérard, B. Legrand, E. Costard, and V. Thierry-Mieg.  Optical study of GaAs/AlAs pillar microcavities with elliptical cross section.  Applied Physics Letters, 72(12):1421–1423, 1998.
  * [29] D. M. Whittaker, P. S. S. Guimaraes, D. Sanvitto, H. Vinck, S. Lam, A. Daraei, J. A. Timpson, A. M. Fox, M. S. Skolnick, Y.-L. D. Ho, J. G. Rarity, M. Hopkinson, and A. Tahraoui.  High Q modes in elliptical microcavity pillars.  Applied Physics Letters, 90(16):161105, 2007.
  * [30] P. Hilaire, C. Antón, C. Kessler, A. Lemaître, I. Sagnes, N. Somaschi, P. Senellart, and L. Lanco.  Accurate measurement of a $96\%$ input coupling into a cavity using polarization tomography.  Appl. Phys. Lett., 112:201101, 2018.
  * [31] P. Androvitsaneas, A. B. Young, J. M. Lennon, C. Schneider, S. Maier, J. J. Hinchliff, G. S. Atkinson, E. Harbord, M. Kamp, S. Höfling, J. G. Rarity, and R. Oulton.  Efficient quantum photonic phase shift in a low Q-factor regime.  ACS Photonics, 6:429–435, 2019.
  * [32] Shuo Sun, Hyochul Kim, Glenn S. Solomon, and Edo Waks.  A quantum phase switch between a single solid-state spin and a photon.  Nature Nanotechnology, 11:539–544, 2016.
  * [33] D.F. Walls and Gerard J. Milburn, editors.  Quantum Optics.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2008.
  * [34] E. Mehdi, M. Gundín, C. Millet, N. Somaschi, A. Lemaître, I. Sagnes, L. Le Gratiet, D. A. Fioretto, N. Belabas, O. Krebs, P. Senellart, and L. Lanco.  Giant optical polarisation rotations induced by a single quantum dot spin.  Nature Communication, 15:598, 2024.
  * [35] Dennis H. Goldstein.  Polarized Light.  CRC Press, 3rd edition, 2017.
  * [36] S. Gerhardt, M. Deppisch, S. Betzold, T. H. Harder, T. C. H. Liew, Ana Predojević, S. Höfling, and C. Schneider.  Polarization-dependent light-matter coupling and highly indistinguishable resonant fluorescence photons from quantum dot-micropillar cavities with elliptical cross section.  Phys. Rev. B, 100:115305, 2019.
  * [37] Maxime Zerbib, Maxime Romanet, Thibaut Sylvestre, Christian Wolff, Birgit Stiller, Jean-Charles Beugnot, and Kien Phan Huy.  Spin-orbit interaction in nanofiber-based brillouin scattering.  Opt. Express, 31(14):22284–22295, 2023.