JFTR is partially supported by the Japan Science and Technology Agency under
CREST Grant Number JPMJCR2014 and by the Japan Society for the Promotion of
Science (JSPS) Grant-in-Aid for Early-Career Scientists under Japan Grant
Number JP23K13012. The authors gratefully acknowledge the anonymous reviewer
for their valuable comments, which significantly enhanced the final article.
Additionally, special thanks to John Sebastian Simon for his correction in the
earlier version of the paper.

## References

  * [ADEK07] L. Afraites, M. Dambrine, K. Eppler, and D. Kateb.  Detecting perfectly insulated obstacles by shape optimization techniques of order two.  Discrete Contin. Dyn. Syst. Ser. B, 8(2), 2007.
  * [ADK08] L. Afraites, M. Dambrine, and D. Kateb.  On second order shape optimization methods for electrical impedance tomography.  SIAM J. Control Optim., 47(3):1556–1590, 2008.
  * [AF03] R. A. Adams and J. J. F. Fournier.  Sobolev Spaces, volume 140 of Pure and Applied Mathematics.  Academic Press, Amsterdam, 2003.
  * [Afr22] L. Afraites.  A new coupled complex boundary method (CCBM) for an inverse obstacle problem.  Discrete Contin. Dyn. Syst. Ser. S, 15(1):23 – 40, 2022.
  * [AMN22] L. Afraites, C. Masnaoui, and M. Nachaoui.  Shape optimization method for an inverse geometric source problem and stability at critical shape.  Discrete Contin. Dyn. Syst. Ser. S, 15(1):1–21, 2022.
  * [AMR09] C. S. Alves, N. F. M. Martins, and N. C. Roberty.  Full identification of acoustic sources with multiple frequencies and boundary measurements.  Inverse Prob. Imaging, 3(2):275–294, 2009.
  * [APR03] G. Alessandrini, L. Del Piero, and L. Rondi.  Stable determination of corrosion by a single electrostatic boundary measurement.  Inverse Problems, 19:973–984, 2003.
  * [AR22] L. Afraites and J. F. T. Rabago.  Shape optimization methods for detecting an unknown boundary with the Robin condition by a single measurement.  Discrete Contin. Dyn. Syst. Ser. S, 2022.
  * [AS07] G. Alessandrini and E. Sincich.  Solving elliptic Cauchy problems and identification of non-linear corrosion.  J. Comput. Appl. Math., 198:307–320, 2007.
  * [Bac09] V. Bacchelli.  Uniqueness for the determination of unknown boundary and impedance with the homogeneous Robin condition.  Inverse Problems, 25:Art. 015004 (4pp), 2009.
  * [BC07] A. Boulkhemair and A. Chakib.  On the uniform Poincaré inequality.  Comm. Partial Differential Equations, 32:1439–1447, 2007.
  * [BNC08] A. Boulkhemair, A. Nachaoui, and A. Chakib.  Uniform trace theorem and application to shape optimization.  Appl. Comput. Math., 7:192–205, 2008.
  * [BNC13] A. Boulkhemair, A. Nachaoui, and A. Chakib.  A shape optimization approach for a class of free boundary problems of Bernoulli type.  Appl. Math., 58:205–221, 2013.
  * [CDK13] F. Caubet, M. Dambrine, and D. Kateb.  Shape optimization methods for the inverse obstacle problem with generalized impedance boundary conditions.  Inverse Problems, 29:Art. 115011 (26pp), 2013.
  * [Che75] D. Chenais.  On the existence of a solution in a domain identification problem.  J. Math. Anal. Appl., 52:189–219, 1975.
  * [CJ99] S. Chaabane and M. Jaoua.  Identification of Robin coefficients by means of boundary measurements.  Inverse Problems, 15:1425–1438, 1999.
  * [CK07] F. Cakoni and R. Kress.  Integral equations for inverse problems in corrosion detection from partial cauchy data.  Inverse Prob. Imaging, 1:229–245, 2007.
  * [CKS10] F. Cakoni, R. Kress, and C. Schuft.  Integral equations for inverse problems in corrosion detection from partial Cauchy data.  Inverse Problems, 26:Art. 095012 24pp., 2010.
  * [Dam02] M. Dambrine.  On variations of the shape Hessian and sufficient conditions for the stability of critical shapes.  Rev. R. Acad. Cienc. Exactas Fis. Nat. Ser. A. Mat., 96:95–121, 2002.
  * [DD97] E. N. Dancer and D. Daners.  Domain perturbation for elliptic equations subject to Robin boundary conditions.  J. Differential Equations, 138:86–132, 1997.
  * [DMNV07] G. Doǧan, P. Morin, R.H. Nochetto, and M. Verani.  Discrete gradient flows for shape optimization and applications.  Comput. Methods Appl. Mech. Engrg., 196:3898–3914, 2007.
  * [DP00] M. Dambrine and M. Pierre.  About stability of equilibrium shapes.  Model Math. Anal. Numer., 34:811–834, 2000.
  * [DSZ03] M. Dambrine, J. Sokolowski, and A. Zochowski.  On analysis in shape optimisation: critical shapes for Neumann problem.  Control Cybern., 30(3):503–528, 2003.
  * [DZ11] M. C. Delfour and J.-P. Zolésio.  Shapes and Geometries: Metrics, Analysis, Differential Calculus, and Optimization, volume 22 of Adv. Des. Control.  SIAM, Philadelphia, 2nd edition, 2011.
  * [EH05] K. Eppler and H. Harbrecht.  A regularized Newton method in electrical impedance tomography using shape hessian information.  Control Cybern., 34:203–225, 2005.
  * [EHS07] K. Eppler, H. Harbrecht, and R. Schneider.  On convergence in elliptic shape optimization.  SIAM J. Control Optim., 46(1):61–83, 2007.
  * [Epp00] K. Eppler.  Boundary integral representations of second derivatives in shape optimization.  Discuss. Math. Differ. Incl. Control. Optim., 20:487–516, 2000\.
  * [Fan22] W. Fang.  Simultaneous recovery of Robin boundary and coefficient for the Laplace equation by shape derivative.  J. Comput. Appl. Math., 413:Art. 114376 13 pp, 2022.
  * [FI07] D. Fasino and G. Inglese.  Recovering nonlinear terms in an inverse boundary value problem for Laplace’s equation: a stability estimate.  J. Comput. Appl. Math., 198:460–470, 2007.
  * [FL04] W. Fang and M. Lu.  A fast collocation method for an inverse boundary value problem,.  Internat. J. Numer. Methods Engrg., 59:1563–1585, 2004.
  * [FLM19] W. Fang, F. Lin, and Y. Ma.  Fast algorithms for boundary integral equations on elliptic domains and related inverse problems.  East Asian J. Appl. Math., 9:485–505, 2019.
  * [FZ09] W. Fang and S. Zeng.  Numerical recovery of Robin boundary from boundary measurements for the Laplace equation.  J. Comput. Appl. Math., 224:573–580, 2009.
  * [GLZ22] Wei Gong, Jiajie Li, and Shengfeng Zhu.  Improved discrete boundary type shape gradients for PDE-constrained shape optimization.  SIAM J. Sci. Comput., 44:A2464–A2505, 2022.
  * [GPT17] M. Giacomini, O. Pantz, and K. Trabelsi.  Certified descent algorithm for shape optimization driven by fully-computable a posteriori error estimators.  ESAIM - Control Optim. Calc. Var., 23:977–1001, 2017.
  * [Gri85] P. Grisvard.  Elliptic Problems in Nonsmooth Domains.  Pitman Publishing, Marshfield, Massachusetts, 1985.
  * [GZ21] Wei Gong and Shengfeng Zhu.  On discrete shape gradients of boundary type for PDE-constrained shape optimization.  SIAM J. Numer. Anal., 59:1510–1541, 2021.
  * [Hec12] F. Hecht.  New development in FreeFem++.  J. Numer. Math., 20:251–265, 2012.
  * [HKKP03] J. Haslinger, T. Kozubek, K. Kunisch, and G. H. Peichl.  Shape optimization and fictitious domain approach for solving free-boundary value problems of Bernoulli type.  Comput. Optim. Appl., 26(3):231–251, 2003.
  * [HKKP04] J. Haslinger, T. Kozubek, K. Kunisch, and G. H. Peichl.  An embedding domain approach for a class of 2-d shape optimization problems: Mathematical analysis.  J. Math. Anal. Appl., 209(2):665–685, 2004.
  * [HM03] J. Haslinger and R. A. E. Mäkinen.  Introduction to Shape Optimization: Theory, Approximation, and Computation.  SIAM, Philadelphia, 2003.
  * [Hol01] L. Holzleitner.  Hausdorff convergence of domains and their boundaries for shape optimal design.  Control Cybern., 30(1):23–44, 2001.
  * [HP18] A. Henrot and M. Pierre.  Shape Variation and Optimization: A Geometrical Analysis, volume 28 of Tracts in Mathematics.  European Mathematical Society, Zürich, 2018.
  * [IK06] O. Ivanyshyn and R. Kress.  Nonlinear integral equations for solving inverse boundary value problems for inclusions and cracks.  J. Integral Equ. Appl., 18(1):13–38, 2006.
  * [IM04] G. Inglese and F. Mariani.  Corrosion detection in conducting boundaries.  Inverse Problems, 20:1207–1215, 2004.
  * [Ing97] G. Inglese.  An inverse problem in corrosion detection.  Inverse Problems, 13:977–994, 1997.
  * [KC98] R. Kress and D. Colton.  Inverse Acoustic and Electromagnetic Scattering Theory, volume 93 of Applied Mathematical Sciences.  Springer, New York, 3rd edition, 1998.
  * [KMI17] T. Kurahashi, K. Maruoka, and T. Iyama.  Numerical shape identification of cavity in three dimensions based on thermal non-destructive testing data.  Eng. Optim., 49(3):434–448, 2017.
  * [KR05] R. Kress and W. Rundell.  Nonlinear integral equations and the iterative solution for an inverse boundary value problem.  Inverse Problems, 21:1207–1223, 2005.
  * [Kre89] Erwin Kreyszig.  Introductory Functional Analysis with Applications.  Wiley, 1989.
  * [KS95] P. G. Kaup and F. Santosa.  Nondestructive evaluation of corrosion damage using electrostatic measurements.  J. Nondestr. Eval., 14:127–136, 1995.
  * [KSV96] P. G. Kaup, F. Santosa, and M. Vogelius.  Method for imaging corrosion damage in thin plates from electrostatic data.  Inverse Problems, 12:279–293, 1996.
  * [Loh87] W. H. Loh.  Modeling and measurement of contact resistances.  Tech. Rep. No. G 830-1, Stanford Electronics Labs., 1987.
  * [Mef17] H. Meftahi.  Stability analysis in the inverse Robin transmission problem.  Math. Meth. Appl. Sci., 40(7):2505–2521, 2017.
  * [MS76] F. Murat and J. Simon.  Sur le contrôle par un domaine géométrique.  Research report 76015, Univ. Pierre et Marie Curie, Paris, 1976.
  * [Neu97] J. W. Neuberger.  Sobolev Gradients and Differential Equations.  Springer-Verlag, Berlin, 1997.
  * [NR00] A. Novruzi and J.-R. Roche.  Newton’s method in shape optimisation: a three-dimensional case.  BIT Numer. Math., 40:102–120, 2000.
  * [Pir84] O. Pironneau.  Optimal Shape Design for Elliptic Systems.  Springer series in Computational Physics. Springer-Verlag, 1984.
  * [PP09] C. D. Pagani and D. Peerotti.  Identifiability problems of defects with Robin condition.  Inverse Problems, 25:Art. 055007 12pp, 2009.
  * [RA18] J. F. T. Rabago and H. Azegami.  Shape optimization approach to defect-shape identification with convective boundary condition via partial boundary measurement.  Japan J. Indust. Appl. Math., 31(1):131–176, 2018.
  * [RA19] J. F. T. Rabago and H. Azegami.  A new energy-gap cost functional cost functional approach for the exterior Bernoulli free boundary problem.  Evol. Equ. Control Theory, 8(4):785–824, 2019.
  * [RA20] J. F. T. Rabago and H. Azegami.  A second-order shape optimization algorithm for solving the exterior Bernoulli free boundary problem using a new boundary cost functional.  Comput. Optim. Appl., 77(1):251–305, 2020.
  * [Run08] W. Rundell.  Recovering an obstacle and its impedance from Cauchy data.  Inverse Problems, 24:1–22, 2008.
  * [Sim80] J. Simon.  Differentiation with respect to the domain in boundary value.  Numer. Funct. Anal. Optim., 2:649–687, 1980.
  * [Sim89] J. Simon.  Second variation for domain optimization problems.  In F. Kappel, K. Kunisch, and W. Schappacher, editors, Control and Estimation of Distributed Parameter Systems, volume 91 of International Series of Numerical Mathematics, pages 361–378, Basel, 1989. Birkhäuser.
  * [Sin10] E. Sincich.  Stability for the determination of unknown boundary and impedance with a Robin boundary condition.  SIAM J. Math. Anal., 42:2922–2943, 2010.
  * [SZ92] J. Sokołowski and J.-P. Zolésio.  Introduction to Shape Optimization: Shape Sensitivity Analysis.  Springer Series in Computational Mathematics. Springer-Verlag, Berlin, Heidelberg, 1992.


# Effective mass suppression in a ferromagnetic two-dimensional electron
liquid

Reza Asgari  School of Physics, Institute for Research in Fundamental
Sciences, (IPM) 19395-5531 Tehran, Iran  T. Gokmen  Department of Electrical
Engineering, Princeton University, Princeton, New Jersey 08544, USA  B.
Tanatar  Department of Physics, Bilkent University, Bilkent, Ankara 06800,
Turkey  Medini Padmanabhan  Department of Electrical Engineering, Princeton
University, Princeton, New Jersey 08544, USA  M. Shayegan  Department of
Electrical Engineering, Princeton University, Princeton, New Jersey 08544, USA

###### Abstract

We present numerical calculations of the electron effective mass in an
interacting, ferromagnetic, two-dimensional electron system. We consider
quantum interaction effects associated with the charge-density fluctuation
induced many-body vertex corrections. Our theory, which is free of adjustable
parameters, reveals that the effective mass is suppressed (relative to its
band value) in the strong coupling limit, in good agreement with recent
experimental results.

###### pacs:

71.10.Ca, 73.20.Mf

##  I Introduction

Two-dimensional electron systems (2DESs) realized at semiconductor interfaces
are of continuing interest ando ; gv_book  from both basic physics and
technological points of view. As a function of the interaction strength, which
is characterized by the ratio $r_{s}$ of the Coulomb energy to Fermi energy,
many novel correlated ground states have been predicted such as a paramagnetic
liquid ($r_{s}<26$), ferromagnetic liquid ($26<r_{s}<35$) and Wigner crystal
($r_{s}>35$) AttaccalitePRL02 . In the paramagnetic liquid phase, interaction
typically leads to an enhancement of effective mass ($m^{*}$) and spin
susceptibility (${\chi}^{*}$ ${\propto}$ $g^{*}m^{*}$), where $g^{*}$ is the
Landé $g^{*}$-factor. Effective mass is an important concept in Landau’s Fermi
liquid theory since it provides a direct measure of the many-body interactions
in the electron system as characterized by increasing $r_{s}$.

The effective mass $m^{*}$ renormalized by interactions has been
experimentally studied SmithPRL72 ; PanPRB99 ; ShashkinPRL03 ; TanPRL05 ;
PadmanabhanPRL08  for various paramagnetic 2DESs as a function of $r_{s}$. In
the highly interacting, dilute, paramagnetic regime ($3<r_{s}<26$), $m^{*}$ is
typically significantly enhanced compared to its band value, $m_{b}$, and
tends to increase with increasing $r_{s}$ SmithPRL72 ; PanPRB99 ;
ShashkinPRL03 ; TanPRL05 ; PadmanabhanPRL08 ; KwonPRB94 ; AsgariSSC04 ;
AsgariPRB05 ; AsgariPRB06 ; GangadharaiahPRL05 ; ZhangPRL05 . A question of
particular interest is the dependence of $m^{*}$ on the 2D electrons’ spin and
valley degrees of freedom as these affect the exchange interaction. Recent
measurements of $m^{*}$ for 2D electrons confined to AlAs quantum wells
revealed that, when the 2DES is fully valley- and spin-polarized, $m^{*}$ is
suppressed down to values near or even slightly below $m_{b}$ PadmanabhanPRL08
; GokmenPRL08 ; GokmenUNP09 . Note that in these experiments, $r_{s}<22$ so
that the 2DES is in the paramagnetic regime, but a strong magnetic field is
applied in order to fully spin-polarize the electrons. Here we present
theoretical calculations indicating that the $m^{*}$ suppression is caused by
the absence (freezing out) of the spin fluctuations. The results of our
$m^{*}$ calculations are indeed in semi-quantitative agreement with the
measurements.

Previous theoretical calculations of the effective mass are mostly performed
within the framework of Landau’s Fermi liquid theory whose key ingredient is
the quasiparticle (QP) concept and its interactions. This entails the
calculation of effective electron-electron interactions which enter the many-
body formalism allowing the calculation of effective mass. A number of works
considered different variants of the leading order in the screened interaction
for the self-energy em ; yarlagadda_1994_2 ; em_bohm ; em_dassarma ; DasPRB04
; zhang ; AsgariSSC04 ; AsgariPRB05 ; AsgariPRB06  from which density, spin-
polarization, and temperature dependence of effective mass are obtained. In
these calculations the on-shell approximation em_bohm ; em_dassarma ; DasPRB04
yields a diverging effective mass but the full solution of the Dyson equation
yields only a mild enhancement. AsgariSSC04 ; AsgariPRB05 ; AsgariPRB06
Almost all these works considered a paramagnetic 2DES as past experiments
concentrated on the effective mass enhancement in partially spin polarized 2D
systems with $r_{s}<26$.

##  II Theory

We consider a ferromagnetic 2DES as a model for a system of electronic
carriers with band mass $m_{b}$ in a semiconductor heterostructure with
dielectric constant $\kappa$. The bare electron-electron interaction is given
by $v_{\bf q}=2\pi e^{2}/(\kappa q)$. At zero temperature there is only one
relevant parameter for the homogeneous, ferromagnetic 2DES, the usual Wigner-
Seitz density parameter $r_{s}=(\pi na^{2}_{B})^{-1/2}$ in which
$a_{B}=\hbar^{2}\kappa/(m_{b}e^{2})$ is the Bohr radius in the medium of
interest.

Figure 1: (color online). Many-body effective mass as a function of $r_{s}$
for $0\leq r_{s}\leq 22$ for a ferromagnetic 2DES.

The QP self-energy with momentum $\bf k$ and frequency $\omega$ in a fully
polarized electron system can be written as

$\Sigma^{\uparrow}({\bf k},\omega)=-\int\frac{d^{2}{\bf q}}{i(2\pi)^{2}}v_{\bf
q}\int_{-\infty}^{\infty}\frac{d\Omega}{2\pi}\frac{1}{\varepsilon({\bf
q},\Omega)}\,\left[\frac{1-n_{\rm F}(\xi^{\uparrow}_{\bf
k})}{\omega+\Omega-\xi^{\uparrow}_{{\bf k}+{\bf q}}/\hbar+i\eta}+\frac{n_{\rm
F}(\xi^{\uparrow}_{\bf k})}{\omega+\Omega-\xi^{\uparrow}_{{\bf k}+{\bf
q}}/\hbar-i\eta}\right]\,.$ (1)

Here $\xi^{\uparrow}_{\bf k}=\varepsilon_{\bf k}-\varepsilon_{F}$ where
$\varepsilon_{\bf k}=\hbar^{2}{\bf k}^{2}/(2m_{b})$ is the single-particle
energy with $\varepsilon_{F}=\hbar^{2}{k^{\uparrow}_{F}}^{2}/(2m_{b})$ and
$k^{\uparrow}_{F}=(4\pi n_{\rm\scriptscriptstyle 2D})^{1/2}=2/(r_{s}a_{B})$,
respectively, being the Fermi energy and wave vector; $n_{\rm F}(k)$ is the
Fermi function. In Eq. (1), $\varepsilon({\bf q},\omega)$ is the dynamical
screening function for which we use the form appropriate for a ferromagnetic
2DES derived from Kukkonen-Overhauser effective interaction. kukkonen  The
many-body exchange and correlation (XC) effects are introduced through the
local-field factors (LFF) $G_{\sigma,\sigma^{\prime}}(q,\omega)$ ($\sigma$ and
$\sigma^{\prime}$ are spin indices) take the Pauli-Coulomb hole around a
charged particle into account. The dynamical screening function reads

$\frac{1}{\varepsilon({\bf q},\omega)}=1+v_{\bf
q}\,\left[1-G^{+}_{\uparrow}({\bf q},\omega)\right]^{2}\,\chi_{\rm\scriptstyle
C}({\bf q},\omega)~{},$ (2)

where $G^{+}_{\uparrow}$ is the LFF associated with charge-fluctuations. This
expression is similar to the Kukkonen and Overhauser interaction kukkonen
where the spin-fluctuation term is dropped. A similar expression has also been
reported in Refs. [ng, ; tanaka, ]. In Eq. (2) $\chi_{\rm\scriptstyle C}({\bf
q},\omega)$ represents the density-density response function, which in turn is
determined by the local-field factor $G^{+}_{\uparrow}({\bf q},\omega)$ via
the relation

$\chi_{\rm\scriptstyle C}({\bf q},\omega)=\frac{\chi^{0}_{\uparrow}({\bf
q},\omega)}{1-v_{\bf q}[1-G^{+}_{\uparrow}({\bf
q},\omega)]\chi^{0}_{\uparrow}({\bf q},\omega)}\,,$ (3)

in which $\chi^{0}_{\uparrow}(q,\omega)$ is the density response function of
the spin-polarized electrons. The expression for the noninteracting density
response function on the imaginary frequency axis is obtained for use in Eq.
(3) as

$\chi^{0}_{\uparrow}({\bf
q},i\Omega)=\frac{m_{b}^{2}}{2\pi\hbar^{2}q^{2}}\left(\sqrt{2}\sqrt{a_{\uparrow}+\sqrt{a_{\uparrow}^{2}+\left(\frac{q^{2}\Omega}{\hbar
m_{b}}\right)^{2}}}-\frac{q^{2}}{m_{b}}\right)\,,$ (4)

where we have defined
$a_{\uparrow}=q^{4}/4m_{b}^{2}-q^{2}{k^{\uparrow}_{F}}^{2}/m_{b}^{2}-\Omega^{2}/\hbar^{2}$.
It is evident that setting $G^{+}_{\uparrow}({\bf q},\omega)=0$, we recover
the standard random phase approximation (RPA). In what follows, we shall make
the common approximation of neglecting the frequency dependence of
$G^{+}_{\uparrow}$.

Quite generally, once the QP retarded self-energy is known, the QP excitation
energy $\delta{\mathcal{E}}^{\uparrow}_{\rm QP}({\bf k})$, which is the QP
energy measured from the chemical potential $\mu^{\uparrow}$ of the
interacting ferromagnetic 2DES, can be calculated by solving self-consistently
the Dyson equation

$\delta{\mathcal{E}}^{\uparrow}_{\rm QP}({\bf k})=\xi^{\uparrow}_{\bf
k}+\left.\Re e\Sigma^{R}_{\rm\scriptstyle ret}({\bf
k},\omega)\right|_{\omega=\,\delta{\mathcal{E}}^{\uparrow}_{\rm QP}({\bf
k})/\hbar}\,.$ (5)

Alternatively, the QP excitation energy can also be calculated from

$\delta{\mathcal{E}}^{\uparrow}_{\rm QP}({\bf k})=\xi^{\uparrow}_{\bf
k}+\left.\Re e\Sigma^{R}_{ret}({\bf
k},\omega)\right|_{\omega=\xi^{\uparrow}_{\bf k}/\hbar}\,.$ (6)

This is called the on-shell approximation (OSA) and it is argued rice  to be a
better approach than solving the full Dyson equation since noninteracting
Green function is used in Eq. (1). Here $\Re e\Sigma^{R}_{ret}({\bf
k},\omega)$ is defined as $\Re e\Sigma^{\uparrow}_{ret}({\bf
k},\omega)-\Sigma^{\uparrow}_{ret}({\bf k^{\uparrow}_{F}},0)$.

The effective mass $m^{*}_{\uparrow}(k)$ is now calculated from

$\frac{1}{m^{*}_{\uparrow}(k)}=\frac{1}{\hbar^{2}k}\frac{d\delta{\mathcal{E}}^{\uparrow}_{\rm
QP}(k)}{dk}\,,$ (7)

where for $\delta{\mathcal{E}}^{\uparrow}_{\rm QP}$ we have at our disposal
the Dyson and OSA approaches. Evaluating $m^{*}_{\uparrow}(k)$ at
$k=k^{\uparrow}_{F}$, one gets the QP effective mass at the Fermi contour.
Clearly from Eqs. (2) and (3) LFF is the basic quantity for an evaluation of
the QP properties. We have used the parameterized forms of LFFs
$G^{+}(q,\zeta)$ and $G^{-}(q,\zeta)$ (and in particular
$G^{+}_{\uparrow}(q)=G^{+}(q,\zeta=1)$ where $\zeta$ is the spin polarization)
of Moreno and Marinescu. moreno

##  III Results and discussion

We now present our numerical results, which are based on the LFF
$G^{+}_{\uparrow}(q)$ as input. In Fig. 1 we show our numerical results of the
QP effective mass both in OSA and Dyson approximations. The QP effective mass
suppression is substantially smaller in the Dyson equation calculation than in
the OSA; the reason is that a significant cancellation occurs between the
numerator and the denominator in the effective mass expression in the Dyson
approach. To clarify the effect of charge-density fluctuation we have also
shown the RPA results which do not take the strong many-body fluctuations into
account. Note that the LFF takes into account multiple scattering events to
infinite order as compared to the RPA where these effects are neglected. In
the limit of small $\zeta$ and $r_{s}\rightarrow 0$, the effective mass can be
analytically shown to be $m^{*}_{\uparrow}/m_{b}=1+(1-\zeta/2.0)r_{s}\ln
r_{s}/(\sqrt{2}\pi)$ which our numerical calculations faithfully reproduce.

Figure 2: (color online). Many-body effective mass as a function of $r_{s}$
for $0\leq r_{s}\leq 22$ for the ferromagnetic 2DES in comparison to
experiments in Ref. [PadmanabhanPRL08, ; GokmenUNP09, ]. Different symbols
denote different samples; triangles: A, squares: B, circles: C, and diamonds:
D.

In Fig. 2 we compare our effective mass calculations with the experimental
results PadmanabhanPRL08 ; GokmenPRL08 ; GokmenUNP09 . The measurements were
made on 2DESs confined to modulation-doped AlAs quantum wells (QWs) of width
4.5, 11, 12, and 15 nm (samples A, B, C, and D). These samples were grown on
GaAs substrates using molecular beam epitaxy. In bulk AlAs, electrons occupy
three degenerate ellipsoidal conduction band valleys at the X-points of the
Brillouin zone with longitudinal and transverse effective masses $m_{l}$=1.05
and $m_{t}$=0.205 (in units of the free electron mass). Thanks to the slightly
larger lattice constant of AlAs compared to GaAs, the AlAs QW layer is under
bi-axial compressive strain. Because of this compression, the 2DES in the
wider QW samples (B, C, and D) occupy two in-plane valleys with their major
axes lying in the plane ShayeganPSS06 . In our measurements on these samples,
we applied uni-axial, in-plane strain to break the symmetry between these two
valleys so that only one in-plane valley, with an $anisotropic$ Fermi contour
and band effective mass of $m_{b}=\sqrt{m_{l}m_{t}}=0.46$ is occupied
ShayeganPSS06 . In sample A, however, thanks to its very small QW width, the
confinement energy of the out-of-plane valley is lower (because of its larger
mass along the growth direction), so that the electrons occupy this valley and
therefore have an $isotropic$ Fermi contour and band effective mass is
$m_{b}=m_{t}=0.205$ ShayeganPSS06 . The effective masses were deduced from the
temperature dependence of the Shubnikov-de Haas oscillations, the details of
which are given in Refs. PadmanabhanPRL08 ; GokmenPRL08 ; GokmenUNP09 . We
emphasize that the data shown here (Fig. 2) were taken on single-valley 2DESs
which were subjected to sufficiently large magnetic fields to fully spin
polarize the electrons .

It appears in Fig. 2 that the OSA accounts overall for the observed reduction
of $m^{*}_{\uparrow}$ below the band value reasonably well. The agreement is
particularly good for the wider samples (B, C, and D) which have $r_{s}>7$.
The $m^{*}$ data for the narrowest sample (A), however, fall above the
theoretical predictions. We do not know the reason for this discrepancy.
However, we point out that, besides the difference in the shapes of the Fermi
contour, there is another difference between sample A and the other three
samples. Because of the very narrow width of sample A’s quantum well and the
prevalence of interface roughness scattering VakiliAPL06 , the mobility of the
electrons in this sample is much lower (about a factor of 6) than in other
samples for comparable $r_{s}$. It is possible that the higher disorder in
sample A is responsible for $m^{*}$ being larger; this conjecture is indeed
consistent with the results of calculations AsgariSSC04  which predict a
larger $m^{*}$ for more disordered samples.

From Figs. 1 and 2 we draw two main conclusions. (i) The RPA and present
results are rather similar in the weak coupling limit ($r_{s}<1$). (ii) In the
strong coupling regime ($r_{s}>3$), however, our theoretical calculations
which incorporate the proper many-body effects exhibit a mass suppression,
similar to the experimental data, while the RPA results show a mass
enhancement and are far from the experimental data. We emphasize that the
effective mass at the Fermi contour is significantly suppressed in the fully
polarized case because of the absence of spin-fluctuation contribution. This
suppression suggests that the anti-symmetric Landau parameter $F^{a}_{1}<0$
and thus higher angular momentum Landau parameters may be negligible in a
fully spin-polarized 2DES.

Figure 3: (color online). Many-body on-shell effective mass as a function of
$k/k_{F}$ at $r_{s}=5$ for 2DES with the combined effect of charge
fluctuations in comparison to paramagnetic 2DES. Figure 4: (color online).
Renormalization constant $Z_{\uparrow}$ as a function of $r_{s}$ for
$0<r_{s}<22$ for a ferromagnetic 2DES.

To gain further insight to the density dependence of $m^{\ast}$, we have
calculated the on-shell effective mass as a function of particle momentum $k$
using Eq. (7) evaluated at $\omega(k)=\xi^{\uparrow}_{\bf k}/\hbar$ and
$r_{s}=5$. More specifically, we use

$\frac{m_{b}}{m^{*}_{{\uparrow}}(k)}=1+\frac{m_{b}}{\hbar^{2}k}\frac{d}{dk}\Re
e\Sigma^{\uparrow}_{\rm ret}(k,\xi^{\uparrow}_{k}),$ (8)

for a ferromagnetic case. The results for both paramagnetic and ferromagnetic
cases are shown in Fig. 3. $m^{*}(k)$ for a paramagnetic 2DES by using
$G^{+}(q,\zeta=0)$ and $G^{-}(q,\zeta=0)$ has a sharp peak around $k\approx
k_{\rm F}$ where $k_{\rm F}=2/(r_{s}a_{B})$ and a resonance like divergent
behavior around $k\approx 2k_{\rm F}$. The peak around $k_{\rm F}$ is
associated with spin fluctuations and the divergent behavior around $2k_{\rm
F}$ is related to density fluctuations. ng ; zhang  In particular, the latter
divergence has been extensively studied by Zhang et al. zhang  within the RPA.
It is related to the dispersion instability and coincides with the plasmon
emission. $m^{*}_{\uparrow}(k)$ for the ferromagnetic 2DES, on the other hand,
clearly shows the disappearance of the peak associated with spin fluctuations.
Thus, $m^{\ast}_{\uparrow}(k)$ is very weakly momentum dependent for
$k<k^{\uparrow}_{F}$, since there is a substantial cancelation between the
residue and the exchange plus line self-energy contribution in this regime
which make the real part of the retarded self-energy approximately linear with
respect to $k$. AsgariPRB05  The divergence associated with charge
fluctuations is still present, showing a negative peak around $k=2k_{\rm F}$.
$m^{\ast}_{\uparrow}(k)$ calculated within the RPA reproduces quantitatively
the divergent behavior associated with charge fluctuations but shows some
structure for $k\leq k^{\uparrow}_{\rm F}$, therefore failing to account for
the absence of spin fluctuations. Our density-dependent effective mass results
(Figs. 1 and 2) are consistent with $m^{\ast}(k)$ calculations which we have
checked for a range of $r_{s}$ values.

We have also calculated the renormalization factor $Z_{\uparrow}(r_{s})$ which
is equal to the discontinuity in the momentum distribution at $k_{F}$ and
defined by $Z^{-1}_{\uparrow}=1-\hbar^{-1}\left.\partial_{\omega}\Re
e\Sigma^{\uparrow}_{\rm ret}(k,\omega)\right|_{k=k^{\uparrow}_{F},\omega=0}$.
The effect of charge-fluctuations is to make the $Z_{\uparrow}$ values larger
at large $r_{s}$ compared to the case when they are not included as shown in
Fig. 4. This means that charge-density fluctuations tend to stabilize the
system, whereas the RPA works in the opposite direction.AsgariPRB05  In the
present case including the LFF helps preserve the Fermi liquid picture in the
low density regime.

We have performed our numerical calculations for strictly 2DES. As indicated
above, experimental samples have a finite thickness in the range of $5-15$ nm.
Our theoretical model may be extended to include the finite quantum well width
effects in the following manner. Choosing, say, an infinite square well model
with width $L$ will modify the bare Coulomb interaction, $v_{\bf q}\rightarrow
v_{\bf q}F(qL)$, in which $F(x)$ is a form factor.gold  For consistency, one
should also calculate the local-field factor $G^{+}_{\uparrow}(q)$ using the
same model for the finite width effects. This would provide a better
comparison with experiments. In our case, the local-field factor we use was
constructedmoreno  by the quantum Monte Carlo data for a strictly 2DES and it
is not straightforward to incorporate the finite with effects within such an
approach. Previous calculationsAsgariPRB06  of the effective mass for a
paramagnetic 2DES suggest that the effect of a finite thickness is to suppress
$m^{*}$. Therefore, we surmise that a similar qualitative effect would occur
for the ferromagnetic 2DES. On the other hand, the finite temperature and
disorder effects have a tendency to enhance the effective massDasPRB04 ;
AsgariSSC04  which may lead to a cancellation. These issues require a more
systematic study.

##  IV Summary

In conclusion, our theoretical calculations incorporating the proper Pauli-
Coulomb hole and multi-scattering processes show that in an interacting, fully
spin-polarized 2DES the absence of spin fluctuations reduces the effective
mass below its band value, in agreement with experimental data. Our results
also demonstrate the inadequacy of RPA to account for the observed effective
mass suppression.

###### Acknowledgements.

R. A. thanks M. Polini for helpful discussions. The work at Princeton
University was supported by the NSF. B. T. is supported by TUBITAK (No.
108T743) and TUBA.

## References

  * (1) T. Ando, A.B. Fowler, and F. Stern, Rev. Mod. Phys. 54, 437-672 (1982) .
  * (2) G.F. Giuliani and G. Vignale, Quantum Theory of the Electron Liquid (Cambridge University Press, Cambridge, England, 2005) .
  * (3) C. Attaccalite, S. Moroni, P. Gori-Giorgi, and G. B. Bachelet, Phys. Rev. Lett. 88, 256601 (2002); B. Tanatar and D.M. Ceperley, Phys. Rev. B 39, 5005 (1989) .
  * (4) J.L. Smith and P.J. Stiles, Phys. Rev. Lett. 29, 102 (1972).
  * (5) W. Pan, D.C. Tsui, and B.L. Draper, Phys. Rev. B 59, 10208 (1999).
  * (6) A.A. Shashkin, S.V. Kravchenko, V.T. Dolgopolov, and T.M. Klapwijk, Phys. Rev. B 66, 073303 (2002).
  * (7) Y.-W. Tan, J. Zhu, H.L. Stormer, L.N. Pfeiffer, K.W. Baldwin, and K.W. West Phys. Rev. Lett. 94, 016405 (2005) .
  * (8) M. Padmanabhan, T. Gokmen, N.C. Bishop, and M. Shayegan, Phys. Rev. Lett 101, 026402 (2008) .
  * (9) Y. Kwon, D.M. Ceperley, and R.M. Martin, Phys. Rev. B 50, 1684 (1994); M. Holzmann, B. Bernu, V. Olevano, R.M. Martin, and D.M. Ceperley, Phys. Rev. B 79, 041308(R) (2009).
  * (10) R. Asgari, B. Davoudi, and B. Tanatar, Solid State Commun. 130, 13 (2004).
  * (11) R. Asgari, B. Davoudi, M. Polini, G.F. Giuliani, M.P. Tosi, and G. Vignale, Phys. Rev. B 71, 045323 (2005) .
  * (12) R. Asgari and B. Tanatar, Phys. Rev. B 74, 075301 (2006) .
  * (13) S. Gangadharaiah and D.L. Maslov, Phys. Rev. Lett. 95, 186801 (2005).
  * (14) Y. Zhang and S. Das Sarma, Phys. Rev. Lett. 95, 256603 (2005).
  * (15) T. Gokmen, M. Padmanabhan, and M. Shayegan, Phys. Rev. Lett 101, 146405 (2008) .
  * (16) T. Gokmen, M. Padhamadnan, K. Vakili, E. Tutuc, and M. Shayegan, Phys. Rev. B 79, 195311 (2009).
  * (17) I.K. Marmorkos and S. Das Sarma, Phys. Rev. B 44, R3451 (1991); H.-J. Schulze, P. Schuck, and N. Van Giai, Phys. Rev. B 61, 8026 (2000) .
  * (18) S. Yarlagadda and G.F. Giuliani, Phys. Rev. B 49, 7887 (1994); 61, 12556 (2000); C.S. Ting, T.K. Lee, and J.J. Quinn, Phys. Rev. Lett. 34, 870 (1975).
  * (19) H. M. Böhm and K. Schörkhuber, J. Phys.: Condens. Matter 12, 2007 (2000) .
  * (20) Y. Zhang and S. Das Sarma, Phys. Rev. B 71, 045322 (2005) .
  * (21) S. Das Sarma, Victor M. Galitski and Ying Zhang, Phys. Rev. B 69, 125334 (2004) .
  * (22) Y. Zhang, V.M. Yakovenko, and S. Das Sarma, Phys. Rev. B 71, 115105 (2005) .
  * (23) C.A. Kukkonen and A.W. Overhauser, Phys. Rev. B 20, 550 (1979) .
  * (24) T.-K. Ng and K.S. Singwi, Phys. Rev. B 34, 7743 (1986) .
  * (25) S. Tanaka and S. Ichimaru, Phys. Rev. B 39, 1036 (1989) .
  * (26) T.M. Rice, Ann. Phys. (N.Y.) 31, 100 (1965) .
  * (27) J. Moreno and D.C. Marinescu, Phys. Rev. B 68, 195210 (2003) .
  * (28) M. Shayegan, E.P. De Poortere, O. Gunawan, Y.P. Shkolnikov, E. Tutuc, K. Vakili, Phys. Stat. Sol. (b) 243, 3629 (2006).
  * (29) K. Vakili, Y.P. Shkolnikov, E. Tutuc, E.P. De Poortere, M. Padmanabhan, and M. Shayegan, Appl. Phys. Lett. 89, 172118 (2006) .
  * (30) A. Gold, Phys. Rev. B 35, 723 (1987).


# Evaporation-induced temperature gradient in a foam column

François Boulogne  Université Paris-Saclay, CNRS, Laboratoire de Physique des
Solides, 91405, Orsay, France.  Emmanuelle Rio  Université Paris-Saclay, CNRS,
Laboratoire de Physique des Solides, 91405, Orsay, France.  Frédéric Restagno
Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405,
Orsay, France.

###### Abstract

Various parameters affect the foam stability: surface and bulk rheology of the
solution, gravitational drainage, mechanical vibrations, bubble gas
composition, and also evaporation. Evaporation is often considered through the
prism of a liquid loss, but it also induces a cooling effect due to the
enthalpy of vaporization. In this study, we combine a theoretical and
experimental approach to explore the temperature field in a foam column
evaporating from the top. We show that a measurable temperature profile exists
in this geometry with temperatures at the interface lower than the
environmental temperature by few degrees. We demonstrate that the temperature
profile is the result of a balance between the enthalpy of vaporization and
heat fluxes originating from the thermal conduction of foam and air and the
thermal radiation. For small foam thicknesses compared to the radius, we found
that the temperature gradient is established over the foam thickness while for
large aspect ratios, the gradient is spanning over a lengthscale comparable to
the tube radius.

##  1 Introduction

Foam stability [1] is one of the crucial parameter for successful applications
[2]. In such an assembly of bubbly, the liquid tend to flow downward due to
gravity (drainage) [3], the small bubbles tend to empty in the bigger ones
(coarsening) [4, 5], neighboring bubbles tend to merge (coalescence) [6], and
liquid evaporation affects the foam wetness. These different mechanisms at the
origin of foam aging are either entangled, coarsening and drainage for
example, or poorly understood like coalescence or evaporation, which has been
mainly ignored in the literature. To build an empirical knowledge, systematic
experimental tests have been developed to compare chemical formulations in
various physical and chemical conditions (pH, temperature, etc.). A common
approach is to use foam columns and measure the temporal evolution of the foam
height [7, 8]. For instance, some commercial instruments rely on the Bikerman
test measuring the equilibrium height in presence of a given bubbling velocity
[9] or on the Ross-Milles test based on the time evolution of the height of a
foam generated beforehand [10].

Most of the time, these tests are performed in closed atmosphere. However, in
ambient conditions, the liquid evaporation plays also a role on the stability
of soap films. This aspect has recently received particular attention on soap
films, based on the consideration that evaporation induces thinning of films
[11, 12, 13]. On foams, the effect of evaporation has also been noticed by
several authors [14, 15, 16, 17]. Li et al. suggested that the film thinning
comes with an increase of the surfactant concentration where evaporation takes
place [18]. Such concentration gradient leads to a Marangoni flow increasing
further the film thinning rate, and thus promoting the film bursting. The
opposite effect, i.e. a stability enhanced by evaporation due to Marangoni
effects has been obtained by Chandran Suja et al. on oil-based systems [19].
As a result, laboratory tests performed with sealed containers may not reflect
the same behavior as in the practical conditions of usage. Therefore, the
effect of evaporation on foams must be carefully considered.

More recently, we have evidenced that the temperature of a soap film is not
always equal to the ambient temperature due to evaporation [20]. The cooling
effect appears to be significant, up to $8$ ∘C, for a film of 12 mm of
diameter and a relative humidity of about 20 %. To the best of our knowledge,
the cooling effect induced by evaporation is overlooked in foams, which may
have important consequences to better understand and predict their stability.

In this paper, we propose to investigate experimentally and theoretically the
cooling induced by the evaporation of a foam column. We measure the
temperature profile in foams of different aspect ratios. Then, we rationalize
these measurements by predicting the steady-state temperature of the foam-air
interface and the resulting temperature profile in the foam column. Finally,
we discuss the predictions offered by the model and we compare the predicted
temperature profiles to our measurements.

##  2 Experimental procedure and observations

Figure 1 is a schematic of the experiment. The column of foam is made in a
cylindrical tube of radius $R=17$ mm and of 20 cm total height. Tubes are in
PMMA (Abaqueplast) with a wall thickness of $3$ mm. A soap solution prepared
by mixing a commercial dish washing soap (Fairy) at 10 %wt with pure water
(resistivity = 18.2 M$\Omega\cdot$cm), is poured in the tube such that the
liquid interface reaches a distance $L_{\rm f}>0$ to the rim. The surface
tension of the liquid is $\gamma=25.4\pm 0.1$ mN/m. The foam is produced by
bubbling air with a pressure generator (OF1, Elveflow, France) through a 32G
blunt needle in the soap solution [21]. This method produces a monodisperse
foam with a bubble diameter of $2b=2$ mm; the bubble size is determined by
image analysis of the rising bubbles [22]. The choice of a monodisperse foam
is motivated by the limited coarsening process [23] and to the seek of a
homogeneous, easy to reproduce material. As illustrated in figure 1, the foam
is generated up to the rim and is placed in an environment controlled in
relative humidity at ${\cal R}_{\rm H}=50$ % [24]. The environmental
temperature is measured for each experiment and is $T_{\infty}=21\pm 1$ ∘C.

Figure 1: Foam column of height $L_{\rm f}$ and radius $R$ evaporating at the
rim, $z=0$, in the atmosphere characterized by the vapor concentration
$c_{\infty}$ at the temperature $T_{\infty}$.

The temperature is measured by a thermocouple probe (type K, NiAl-NiCr,
diameter 0.2 mm, RS PRO) connected to a digital thermometer (RS PRO 1314). The
probe is directly plunged in the foam. The position of the probe is controlled
by a motorized translation stage (LTS150, Thorlabs) and the temperature is
recorded once stabilized at each position. Temperatures are measured with a
typical uncertainty of $\pm 0.1$ ∘C and the positions with a typical
uncertainty of 1 mm. The relative humidity is measured with a HIH-4021-003
sensor from Honeywell, USA.

Typical measurements are presented in figure 2 where we calculate the
difference of temperature between the thermocouple at a given position $T(z)$
and the ambient temperature $T_{\infty}$ measured with the same instrument. We
show two temperature profiles for foam thicknesses smaller and larger than the
tube radius, i.e. $L_{\rm f}/R=0.6$ and $L_{\rm f}/R=3.9$ respectively. Within
the duration of the experiments (typically 10 minutes), we have not noticed a
significant foam aging or bubble bursting that would have changed the typical
thickness of the foam layer or the bubble size.

The temperatures of the foam/atmosphere interface are nearly equal in both
situations, to the resolution of the temperature measurements and the
determination of the position of the interface ($\pm$ 1 mm). This surface
temperature is $2.5$ ∘C lower than the environmental temperature for both foam
heights. However, the characteristic lengthscale of the temperature variation
differs with the foam height. For an aspect ratio $L_{\rm f}/R\simeq 0.6$, the
temperature nearly reaches the ambient temperature at the foam-liquid
interface, which results in a sharp temperature increases with the penetration
in the foam. In contrast, for the thicker foam with an aspect ratio $L_{\rm
f}/R\simeq 3.9$, the temperature variation is smoother and the ambient
temperature is reached ahead of the liquid-foam interface.

Additionally, we checked by weighing samples in identical conditions that the
soap solution evaporates at the same rate as pure water, indicating that the
solutes do not alter significantly the chemical potential of the solution
[20].

Figure 2: Measurement of the temperature profile for two foam thicknesses at a
relative humidity ${\cal R}_{\rm H}=48$ %. These thicknesses are materialized
by vertical dashed lines at $-z=L_{\rm f}=11$ and 66 mm. The column radius is
$17$ mm, such that the ratios $L_{\rm f}/R$ are 0.6 and 3.9, respectively.

In the next section, we aim at predicting the temperature profile in the foam,
which includes the determination of the temperature of the interface due to
evaporation.

##  3 Model

###  3.1 Evaporative flux

The atmosphere is characterized by its temperature $T_{\infty}$ and it
relative humidity $\mathcal{R}_{\rm H}$, which is defined as
$p_{\infty}/p_{\rm sat}(T_{\infty})$ where $p_{\infty}$ is the partial vapor
pressure at $T_{\infty}$ and $p_{\rm sat}(T_{\infty})$, the saturated vapor
pressure of water. The temperature variation of the saturated vapor pressure
is well described by the Antoine equation $p_{\rm
sat}(T)=p^{\circ}\,10^{A-B/(C+T)}$ where $p^{\circ}=10^{5}$ Pa and the
coefficients where $A=5.341\pm 0.003$ K, $B=1807.5\pm 1.6$ K, and $C=-33.9\pm
0.1$ K are fitted with data from [25].

In a steady state regime, assuming diffusive transfers between the foam and
the atmosphere, the vapor concentration field is the solution of the Laplace
equation. The total evaporative flux $Q_{\rm ev}$ of a circular disk of radius
$R$ is proportional to the difference of vapor concentrations $\Delta
c^{\star}$ between the environment at a temperature $T_{\infty}$ and above the
interface at a temperature $T_{\rm i}$. This evaporative flux writes

$Q_{\rm ev}=4{\cal D}R\Delta c^{\star},$ (1)

where ${\cal D}$ is the diffusion coefficient of vapor in air [26, 27].

The difference of vapor concentrations $\Delta c^{\star}=c_{\infty}-c_{\rm
sat}(T_{\rm i})$ can be related to the difference of vapor pressures. Denoting
$P$ the atmospheric pressure, $\rho_{\rm air}$ the air density, $M_{\ell}$ and
$M_{\rm air}$ the molar weight of the liquid and air respectively, we have
$\Delta c^{\star}\simeq\frac{\rho_{\rm air}M_{\ell}}{M_{\rm air}}\frac{\Delta
p^{\star}}{P}$ with $\Delta p^{\star}=p_{\infty}-p_{\rm sat}(T_{\rm i})$.

Due to the enthalpy of vaporization, the temperature of the interface $T_{\rm
i}$ is lower than the ambient temperature. Thus, the interface receives a heat
flux from the environment. Three contributions can be identified: by
conduction from the surrounding atmosphere, by conduction from the foam, and
also by radiation. In the following paragraphs, we model these heat transfers.

###  3.2 Heat fluxes

#### Heat flux from the atmosphere

Similarly to the vapor concentration field, the temperature field is the
solution of a Laplace equation with $T_{\rm i}$ at the interface and
$T_{\infty}$ far from it. Thus, the heat flux from the atmosphere is

$Q_{\rm h1}=4\lambda_{\rm air}R\Delta T^{\star},$ (2)

where $\lambda_{\rm air}=0.025$ W$\cdot$K${}^{-1}\cdot$m-1 is the thermal
conductivity of air and $\Delta T^{\star}=T_{\infty}-T_{\rm i}$ is the
temperature difference. Next, to determine the heat flux from the foam, we
must comment its thermal conductivity.

#### Thermal conductivity of foams

The foam is supposed to be a continuous medium, which is valid if the
lengthscale of the heterogeneities is much smaller than the size of the
material, i.e. for $b\ll\\{R,L_{\rm f}\\}$. The foam and the soap solution are
characterized by their thermal conductivity, $\lambda_{\rm f}$ and
$\lambda_{\ell}$, respectively. Leach proposed a model combining series and
parallel contributions to the thermal conductivity [28]. At the first leading
order in $\varphi$, the liquid fraction, the thermal conductivity of foam can
be estimated as

$\lambda_{\rm f}=\lambda_{\rm air}+\frac{2}{3}\lambda_{\ell}\varphi,$ (3)

with $\lambda_{\ell}=0.61$ W$\cdot$K${}^{-1}\cdot$m-1.

The liquid fraction of the foam results from the balance between gravity and
capillary suction, which implies that the liquid fraction decreases with the
altitude. The liquid fraction profile $\varphi(z)$ is obtained by solving the
drainage equation in a stationary regime (See Eq. (3.89) in [29]), which gives

$\varphi(z)=\hat{\varphi}\left(\frac{L_{\rm
f}+z}{\ell_{c}}+\left(\frac{\varphi_{c}}{\hat{\varphi}}\right)^{-1/2}\right)^{-2},$
(4)

where $\varphi_{c}=0.26$ is the fraction of gaps in a close-packing of hard
spheres and $\hat{\varphi}=\ell_{c}^{2}/b^{2}\delta^{2}$ with the capillary
length $\ell_{c}=\sqrt{\gamma/\rho g}$ and a geometric constant $\delta=1.73$.

Figure 3: The main plot shows the thermal conductivity of the foam as the
function of the altitude, computed from equation 3. As an indication, the
thermal conductivity values of water and air are also represented. In the
inset, the liquid fraction is plotted as a function of the vertical position
from equation 4. Computations are performed for a bubble radius $b=1$ mm and
$\ell_{c}=1.6$ mm. As an indication, at $L_{\rm f}+z=25$ mm, $\lambda_{\rm
f}/\lambda_{\rm air}=1.04$.

Equation 4 is plotted in the inset of figure 3 for a bubble radius $b=1$ mm. A
wet region located near the liquid-foam interface spans over a few bubble
layers where a sharp decrease of the liquid fraction occurs, followed by a
smoother decrease with the altitude [29, 30]. From the liquid fraction
profile, we plot in figure 3 the thermal conductivity as a function of the
altitude $z$ with equation 3. Again, we notice a sharp variation of
$\lambda_{\rm f}$ over the first few layers of bubbles and the thermal
conductivity converges quickly to the air conductivity.

In the following, we consider that the thermal conductivity of the foam is
independent of the vertical position $z$ and we take $\lambda_{\rm
f}\simeq\lambda_{\rm air}$. This approximation allows to perform analytical
calculations that will facilitate the interpretation of the predictions made
by the model. After a comparison of the predictions with experimental
measurements, we will comment further this approximation in the discussion.

#### Heat flux from the foam column

We propose to describe the temperature field in the center of the foam $T(z)$.
The environment is at a temperature $T_{\infty}$. Due to the small thickness
of the tube wall compared to the tube radius and the larger thermal
conductivity of the wall compared to the foam, the temperature variation in
the tube wall is neglected. Thus, the temperature profile $T(z)$ at the center
of the column results of the heat transfers from the column periphery, at
$T_{\infty}$, the interface at $T_{\rm i}$, and the liquid, which is far from
the liquid-foam interface at $T(z\rightarrow-\infty)=T_{\infty}$.

A heat flux balance on a slice between $z$ and $z+{\rm d}z$ yields the
differential equation

$\frac{{\rm d}^{2}T}{{\rm
d}z^{2}}+\frac{2}{R^{2}}\left(T_{\infty}-T(z)\right)=0.$ (5)

This equation is valid both in the foam, which lies from $z=0$ to $z=-L_{\rm
f}$, and in the liquid from $z=-L_{\rm f}$ to $z\rightarrow-\infty$. The
boundary conditions associated with the differential equation 5 are

$\displaystyle T(0)$ $\displaystyle=T_{\rm i},$ (6a) $\displaystyle
T(z\rightarrow-\infty)$ $\displaystyle=T_{\infty},$ (6b) $\displaystyle
T\left(-L_{\rm f}^{-}\right)$ $\displaystyle=T\left(-L_{\rm f}^{+}\right),$
(6c) $\displaystyle\lambda_{\rm f}\left.\frac{{\rm d}T}{{\rm
d}z}\right|_{-L_{\rm f}^{+}}$ $\displaystyle=\lambda_{\rm\ell}\left.\frac{{\rm
d}T}{{\rm d}z}\right|_{-L_{\rm f}^{-}},$ (6d)

where equations 6a and 6b are the temperature conditions at the interface and
in the liquid far from the foam, respectively; equation 6c corresponds to the
temperature continuity at the liquid-foam interface, and 6d is the continuity
of the thermal flux at the liquid-foam interface.

The solution of equation 5 is expressed in the foam and the liquid
respectively, and reads

$T_{\ell,\rm f}(z)=T_{\infty}-\Delta T^{\star}\left[\alpha_{\ell,\rm
f}e^{-\sqrt{2}z/R}-\beta_{\ell,\rm f}e^{\sqrt{2}z/R}\right],$ (7)

where the coefficients $\alpha_{\rm f}$, $\beta_{\rm f}$, $\alpha_{\ell}$,
$\beta_{\ell}$ are determined with the boundary conditions 6a-d. Let us denote
$k^{\pm}=\exp(\pm\sqrt{2}L_{\rm f}/R)$ and $\Lambda=\lambda_{\rm
f}/\lambda_{\ell}$. Then, we have

$\displaystyle\alpha_{\ell}$ $\displaystyle=0,$ (8a)
$\displaystyle\beta_{\ell}$ $\displaystyle=\frac{k^{+}}{k^{-}}+\beta_{\rm
f}\left(1+\frac{k^{+}}{k^{-}}\right),$ (8b) $\displaystyle\alpha_{\rm f}$
$\displaystyle=1+\beta_{\rm f},$ (8c) $\displaystyle\beta_{\rm f}$
$\displaystyle=\frac{k^{+}(1+\Lambda)}{k^{-}(1-\Lambda)-k^{+}(1+\Lambda)}.$
(8d)

The heat flux from the bulk of the foam to the interface is $Q_{\rm h2}=-\pi
R^{2}\lambda_{\rm
f}\left.\overrightarrow{z}\cdot\overrightarrow{\nabla}T\right|_{z=0}$.
Combined with equation 7, we obtain

$Q_{\rm h2}=-\pi\sqrt{2}R\lambda_{\rm f}\left(1+2\beta_{\rm f}\right)\Delta
T^{\star}.$ (9)

It is convenient to define the total heat flux by conduction $Q_{\rm h}=Q_{\rm
h1}+Q_{\rm h2}$, which can be written in the form

$Q_{\rm h}=\lambda_{\rm eff}R\Delta T^{\star},$ (10)

with an effective thermal conductivity $\lambda_{\rm eff}=4\lambda_{\rm
air}-\pi\sqrt{2}\left(1+2\beta_{\rm f}\right)\lambda_{\rm f}$. This effective
thermal conductivity depends on foam aspect ratio $L_{\rm f}/R$ through the
coefficient $\beta_{\rm f}$ as shown in figure 4. In the limit of large aspect
ratios, we have $\beta_{\rm f}\rightarrow-1$. Thus, the effective thermal
conductivity becomes independent of the aspect ratio and the value is
$\lambda_{\rm eff}^{\rm lim}=4\lambda_{\rm air}+\pi\sqrt{2}\lambda_{\rm f}$.
Both terms in the expression of $\lambda_{\rm eff}$ are of the same order of
magnitude, which reflects the significance of both $Q_{\rm h1}$ and $Q_{\rm
h2}$ in the heat exchange with the interface.

Figure 4: The effective thermal conductivity $\lambda_{\rm eff}$ as defined in
equation 10 is plotted as a function of the foam aspect ratio $L_{\rm f}/R$ in
orange. The dashed black line is the limit for large aspect ratio
$\lambda_{\rm eff}^{\rm lim}=4\lambda_{\rm air}+\pi\sqrt{2}\lambda_{\rm f}$.

#### Radiative flux

Radiation is known to be significant to estimate the cooling effect when the
characteristic size of the evaporating surface is typically larger than
several millimeters [31, 20]. We describe the radiative flux bringing heat to
the surface by the Stefan-Boltzmann equation

$Q_{\rm rad}=\pi R^{2}\epsilon\sigma(T_{\infty}^{4}-T_{\rm i}^{4}),$ (11)

where $\epsilon$ is the emissivity and $\sigma$ the Stefan-Boltzmann constant.
Below, we will comment the significance of this radiative flux with respect to
the conductive heat flux $Q_{\rm h}$.

Now that we determined the temperature profile in the column and the heat
transfers, the temperature of the interface $T_{\rm i}$ remains to be
calculated to close the problem. In the next paragraph, we derive this
temperature from the balance between the energy absorbed by the evaporation
and the heat transfers.

###  3.3 Energy balance

In a steady state regime, the energy balance writes

$h_{\rm ev}Q_{\rm ev}=-Q_{h}\left(1+\frac{Q_{\rm rad}}{Q_{h}}\right),$ (12)

where $h_{\rm ev}$ is the enthalpy of vaporization. For small temperature
differences, i.e. $|T_{\rm i}-T_{\infty}|/T_{\infty}\ll 1$, the ratio of
radiative and conductive fluxes can be simplified as

$\frac{Q_{\rm rad}}{Q_{\rm h}}=\frac{\pi R\epsilon\sigma(T_{\infty}^{4}-T_{\rm
i}^{4})}{\lambda_{\rm eff}(T_{\infty}-T_{\rm i})}\simeq\frac{4\pi
R\epsilon\sigma T_{\infty}^{3}}{\lambda_{\rm eff}}.$ (13)

Considering large foam aspect ratios, the radiative flux appears to be as
significant as the conductive flux for a critical radius $R_{\rm
c}=\lambda_{\rm eff}^{\rm lim}/(4\pi\epsilon\sigma T_{\infty}^{3})$. At
$21~{}^{\text{o}}$C, we find $R_{\rm c}\simeq 7$ mm, which indicates that
$Q_{\rm rad}$ cannot be neglected in our experiments.

From equation 12, with equations 1, 10, and 13, we relate the difference of
vapor pressures $\Delta p^{\star}$ to the temperature difference $\Delta
T^{\star}$ that reads

$\Delta p^{\star}=-P{\cal A}\Delta T^{\star},$ (14)

where ${\cal A}$ is the so-called psychrometer coefficient, which is for this
system,

${\cal A}=\frac{M_{\rm air}}{\rho_{\rm air}M_{\ell}}\frac{\lambda_{\rm
eff}}{4h_{\rm ev}{\cal D}}\left(1+\frac{4\pi R\epsilon\sigma
T_{\infty}^{3}}{\lambda_{\rm eff}}\right).$ (15)

We can remark that $\Delta p^{\star}$ depends on the temperature of the
interface $T_{\rm i}$, which is necessary to take into account the variation
of the evaporative flux with the temperature of the interface. The temperature
dependence of $\Delta p^{\star}$ being non-linear due to Antoine equation, the
interfacial temperature $T_{\rm i}$ cannot be determined analytically from
equation 14.

As a consequence, for a given column radius $R$, ambient temperature
$T_{\infty}$, and relative humidity ${\cal R}_{\rm H}=(p_{\rm sat}(T_{\rm
i})-\Delta p^{\star})/p_{\rm sat}(T_{\infty})$, we seek numerically for the
temperature $T_{\rm i}$ from equation 14 with the Antoine equation by using a
Newton procedure [32]. Next, the temperature profile is computed with
equations 7 and 8.

From equation 15, it is worth noting that the psychrometer coefficient depends
on lengthscales through (a) $\lambda_{\rm eff}$, which depends itself of
$L_{\rm f}/R$ for small aspect ratios, and (b) directly of the tube radius $R$
due to the ratio $Q_{\rm rad}/Q_{\rm h}$ (Eq. 13). Thus, the psychrometer
coefficient is a function of $(R,L_{\rm f}/R)$. Therefore, the temperature of
the interface cannot be rescaled neither by the radius $R$, nor by the aspect
ratio $L_{\rm f}/R$, and so the temperature profiles.

In the next section, we discuss the results of this model, and we compare the
predictions to experiments.

##  4 Discussion

Before comparing the predictions to the experiments, we start by analyzing the
influence of the model parameters on the temperature of the interface. In
figure 5(a), we plot the cooling effect as a function of the foam column
radius for different values of the relative humidity and a constant foam
aspect ratio.

Evaporation and thus the cooling effect being driven by the relative humidity,
a dry atmosphere leads to a lower temperature of the interface. This statement
is a direct consequence of the psychrometric equation (Eq. 14). In addition,
the increase of the foam column radius reduces the cooling effect. Indeed, in
figure 5(a), the ratio $Q_{\rm rad}/Q_{\rm h}$ increases linearly with the
radius (Eq. 13) such that wider columns receive a more important contribution
of the radiative transfer, which hinders the cooling effect.

The foam aspect ratio also has an effect on the temperature of the interface.
In figure 5(b), we plot the temperature difference as a function of the radius
for aspect ratio smaller and larger than unity. For small aspect ratios,
typically $L_{\rm f}/R<1$, the temperature of the interface decreases with the
foam thickness at a given radius. This observation is directly linked to the
increased effective thermal conductivity $\lambda_{\rm eff}$ observed in
figure 4. Schematically, a small foam thickness reduces the insulation of the
foam-vapor interface with the underlying liquid, which has a significantly
larger thermal conductivity than air. Thus, increasing the aspect ratio beyond
unity has a little effect on the interface temperature, as the temperature
profile becomes dominated by the heat flux from the wall of the container.

Figure 5: Theoretical prediction of the cooling effect on the interface as a
function of the foam column radius. (a) Effect of the relative humidity for
$L_{\rm f}/R=3$. (b) Influence of the foam height with respect to the foam
radius. The ambient temperature is set at $T_{\infty}=20$ ∘C and ${\cal
R}_{\rm H}=50$ %.  Figure 6: Comparison between the measurement of the
temperature profile for various foam thicknesses with the prediction given by
equation 7. The experiments are performed in a column of radius $R=17$ mm at a
relative humidity ${\cal R}_{\rm H}=0.48$. The vertical dashed lines
correspond to the positions of the liquid/foam interface $L_{f}$.

Now that the effects of the relative humidity and the size of the tube are
clarified, we can focus on the temperature profile in the foam for different
aspect ratios. In figure 6, we show the experimental measurements presented in
figure 2 with additional foam heights. We also plot the temperature profiles
given by equation 7 without any fitting parameter. We observe that the model
is in excellent agreement with the measurements, within the experimental
uncertainties.

The temperature profiles presented in figure 6 clearly illustrate the
insulating properties of foam, especially for small aspect ratios where a neat
rupture in the profile is observed at the liquid-foam interface. In this case,
the lengthscale associated with the temperature gradient is the foam
thickness. Also, for larger aspect ratios, the profiles tend to collapse on a
single curve. Again, this collapse indicates that the temperature is mainly
the results of the radial heat flux over the vertical flux. Thus, for large
foam aspect ratios, the characteristic lengthscale of the temperature gradient
is the tube radius.

Finally, we come back to the approximation we made on the thermal conductivity
of the liquid foam. The assumption of a uniform thermal conductivity to
describe the temperature profile in the foam columns is well validated by the
experiments as shown in figure 6. Indeed, we expect that this approximation is
less satisfactory when the temperature gradient, and thus the heat flux, at
the liquid-foam interface becomes significant, which is the case for thin foam
layers. Temperatures of thinner foam layers must deviate from the present
model, but this is difficult to render with our experimental protocol as the
foam thickness will be composed of only a few bubble layers, conflicting with
our continuum approach.

##  5 Conclusions

In this work, we have investigated experimentally and theoretically the
temperature profile that is established in a foam column evaporating from the
top. Experimentally, we observed that a temperature profile is perfectly
measurable with temperature differences between the interface and the ambient
temperature of few degrees. Considering the foam aspect ratio as the ratio
between the foam thickness and the tube radius, we found for thin foams, the
lengthscale of the temperature gradient is set by the foam thickness.
Conversely, for large aspect ratios, the lengthscale is the tube radius.

We successfully modeled the temperature profile in the column with an
analytical expression by considering, in a 1D approximation, a heat flux
balance between the radial and the axial directions, assuming a constant
thermal conductivity in the foam layer. The temperature of the interface
results from an energy balance between the energy associated with the enthalpy
of vaporization and the thermal flux from the environment. The latter
originates from three contributions: the heat conduction by the foam and by
the surrounding air, as well as the radiative flux. The good comparison of our
theoretical predictions with our experiments validates the aforementioned
approximations.

Further studies will be necessary to evaluate how the temperature gradient
combined with the evaporation plays a role on the foam stability. The foam
drying due to evaporation promotes bubble bursting [33], while the role of
thermal gradients is more subtle. For instance, evaporation can increase the
concentration of non-volatile compounds, increasing also the bubble stability
[34]. a Marangoni flow in surface bubbles placed in a cool atmosphere has been
found to enhance the stability by thickening the soap film [35], which remains
an open question on foams.

## Acknowledgments

The authors thank A. Commereuc and M. Pasquet for stimulating discussions. The
authors acknowledge for funding support from the French Agence Nationale de la
Recherche in the framework of project AsperFoam - 19-CE30-0002-01.

## References

  * [1] D. Langevin.  Foam stabilization mechanisms.  In Foam Films and Foams, pages 213–232. CRC Press, 2018.
  * [2] C. Hill and J. Eastoe.  Foams: From nature to industry.  Adv. Colloid Interface Sci., 247:496–513, 2017.
  * [3] S. Cohen-Addad, R. Höhler, and O. Pitois.  Flow in foams and flowing foams.  Annu. Rev. Fluid Mech., 45:241–267, 2013.
  * [4] A. Saint-Jalmes.  Physical chemistry in foam drainage and coarsening.  Soft Matter, 2(10):836–849, 2006.
  * [5] P. Stevenson.  Inter-bubble gas diffusion in liquid foam.  Current Opinion in Colloid & Interface Science, 15(5):374–381, 2010.
  * [6] E. Rio and A.-L. Biance.  Thermodynamic and mechanical timescales involved in foam film rupture and liquid foam coalescence.  Chem. Phys. Chem., 15(17):3692–3707, 2014.
  * [7] R. J. Pugh.  Bubble and foam chemistry.  Cambridge University Press, 2016.
  * [8] J. Wang, A. V. Nguyen, and S. Farrokhpay.  A critical review of the growth, drainage and collapse of foams.  Adv. Colloid Interface Sci., 228:55–70, 2016.
  * [9] J.J. Bikerman.  The unit of foaminess.  Trans. Faraday Soc., 34:634–638, 1938.
  * [10] J. Ross and G. D. Miles.  An apparatus for comparison of foaming properties of soaps and detergents.  Oil & Soap, 18(5):99–102, 1941.
  * [11] L. Champougny, J. Miguet, R. Henaff, F. Restagno, F. Boulogne, and E. Rio.  Influence of Evaporation on Soap Film Rupture.  Langmuir, 34(10):3221–3227, 2018.
  * [12] S. Poulain, E. Villermaux, and L. Bourouiba.  Ageing and burst of surface bubbles.  J. Fluid Mech., 851:636–671, 2018.
  * [13] M. Pasquet, F. Boulogne, J. Saint-Anna, F. Restagno, and E. Rio.  Impact of physical-chemistry on the film thinning in surface bubbles.  Soft Matter, 18:4536–4542, 2022.
  * [14] A. V. Pertsov, V. F. Borachuk, B’ E. Chlstyakov, and E. D. Shchukin.  Evaporation of the dispersion medium and foam stability.  Dokl. Akad. Nauk SSSR, 1980.
  * [15] D. R. Ekserova and P. M. Krugliakov.  Foam and foam films: theory, experiment, application.  Number vol. 5 in Studies in interface science. Elsevier, Amsterdam ; New York, 1998.
  * [16] P. Stevenson.  The wetness of a rising foam.  Ind. Eng. Chem. Res., 45(2):803–807, January 2006.
  * [17] S. T. Tobin, A. J. Meagher, B. Bulfin, M. Möbius, and S. Hutzler.  A public study of the lifetime distribution of soap films.  Am. J. Phys., 79(8):819–824, August 2011.
  * [18] X. Li, S. I. Karakashev, G. M. Evans, and P. Stevenson.  Effect of environmental humidity on static foam stability.  Langmuir, 28(9):4060–4068, 2012.
  * [19] V. Chandran Suja, A. Kar, W. Cates, S.M. Remmert, P.D. Savage, and G.G. Fuller.  Evaporation-induced foam stabilization in lubricating oils.  Proc. Natl. Acad. Sci., 115(31):7919–7924, 2018.
  * [20] F. Boulogne, F. Restagno, and E. Rio.  Measurement of the temperature decrease in evaporating soap films.  Phys. Rev. Lett., 129:268001, 2022.
  * [21] W. Drenckhan and A. Saint-Jalmes.  The science of foaming.  Adv. Colloid Interface Sci., 222:228 – 259, 2015.
  * [22] M. Marchand, F. Restagno, E. Rio, and F. Boulogne.  Roughness-induced friction in liquid foams.  Phys. Rev. Lett., 124:118003, 2020.
  * [23] W. Drenckhan and D. Langevin.  Monodisperse foams in one to three dimensions.  Curr. Opin. Colloid Interface Sci., 15(5):341–358, 2010.
  * [24] F. Boulogne.  Cheap and versatile humidity regulator for environmentally controlled experiments.  Eur. Phys. J. E, 42(4):51, 2019.
  * [25] D.W.H. Rankin.  Handbook of chemistry and physics.  Taylor & Francis, 2009.
  * [26] J. R. Cooke.  Some theoretical considerations in stomatal diffusion: A field theory approach.  Acta Biotheor., 17(3):95–124, 1967.
  * [27] N. N. Lebedev.  Special functions and their applications.  Courier Corporation, 1965.
  * [28] A.G. Leach.  The thermal conductivity of foams. I. Models for heat conduction.  J. Phys. D: Appl. Phys., 26(5):733–739, may 1993.
  * [29] I. Cantat, S. Cohen-Addad, F. Elias, F. Graner, R. Höhler, O. Pitois, F. Rouyer, A. Saint-Jalmes, R. Flatman, and S. Cox.  Foams: Structure and Dynamics.  OUP Oxford, 2013.
  * [30] A. Maestro, W. Drenckhan, E. Rio, and R. Hohler.  Liquid dispersions under gravity: volume fraction profile and osmotic pressure.  Soft Matter, 9:2531–2540, 2013.
  * [31] L. E. Hill, O. W. Griffith, and M. Flack.  V. The measurement of the rate of heat-loss at body temperature by convection, radiation, and evaporation.  Philosophical Transactions of the Royal Society of London. Series B, Containing Papers of a Biological Character, 207(335-347):183–220, 1916.
  * [32] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, I. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [33] V. Carrier and A. Colin.  Coalescence in draining foams.  Langmuir, 19(11):4535–4538, 2003.
  * [34] A. Roux, A. Duchesne, and M. Baudoin.  Everlasting bubbles and liquid films resisting drainage, evaporation, and nuclei-induced bursting.  Phys. Rev. Fluids, 7:L011601, Jan 2022.
  * [35] S. Nath, G. Ricard, P. Jin, A. Bouillant, and D. Quéré.  Thermal marangoni bubbles.  Soft Matter, 18:7422–7426, 2022.

See pages 1- of SI.pdf