# Origin of the Apparent Electric-Field Dependence of Electrostrictive
Coefficients

###### Abstract

Electrostrictive materials exhibit a strain that is proportional to the square
of the induced polarization. In linear dielectrics where the permittivity is
constant, this electromechanical strain is also proportional to the square of
the electric field. However, under increasing amplitudes of the driving field,
the electromechanical strain sometimes saturates; the electrostrictive
coefficients therefore appear to depend on the amplitude of the electric field
used to measure them. Here, we present a methodology showing that this
apparent field dependence is a consequence of neglecting higher-order
electromechanical phenomena. When these are taken into account, not only do
the electrostrictive coefficients remain constant but the signs of the high-
order coefficients enable the prediction of the saturation behavior from a
single measurement. We illustrate this approach on both classical and non-
classical (so-called “giant”) electrostrictors.

###### keywords:

Electrostriction, Dielectrics, Electromechanical Materials

Jiacheng Yu Abdelali Zaki Killian Mache Omar Ibder Sandrine Coste Maud Barré
Philippe Lacorre Pierre-Eymeric Janolin*

Dr. Jiacheng Yu, Dr. Abdelali Zaki, Omar Ibder, Prof. Pierre-Eymeric Janolin*
Université Paris-Saclay, CNRS, CentraleSupélec, laboratoire SPMS, 8-10 rue
Joliot Curie, 91190 Gif-sur-Yvette, France
Université Paris-Saclay, CNRS, CentraleSupélec, SPMS-Pytheas Technology SPyMS
LabCom, 8-10 rue Joliot Curie, 91190 Gif-sur-Yvette, France
Email Address<EMAIL_ADDRESS>
Killian Mache, Dr. Sandrine Coste, Dr. Maud Barré, Dr. Philippe Lacorre
Institut des Molécules et Matériaux du Mans, IMMM-UMR CNRS 6283, Le Mans
Université, Avenue Olivier Messiaen, F-72085 Le Mans Cedex 9, France

##  1 Introduction

Within the last decade, the interest in electrostrictive materials has gained
momentum after the discovery of “giant” or “non-conventional” electrostrictors
[1, 2, 3, 4, 5], exhibiting electromechanical coupling constants at least ten
times larger than expected from the empirical proxy[6]. Recently, giant
effective piezoelectric constants (up to 200 000 pm/V at 0.1 Hz) have been
reported in thin films of Gd-doped ceria[7]. Electrostriction is of interest
for electromechanical applications both as actuators and sensors. As it does
not rely on the existence of ferroelectric domains, it is expected to exhibit
lower losses and superior fatigue performances. A better understanding and
quantification of the fundamental mechanism of the electrostrictive response
will enable improvements to strengthen its potential for applications.
Electrostriction belongs to the class of nonlinear electromechanical couplings
that develop under intense electric fields. Such fields often correspond to
moderate voltages applied to thin films or to larger samples under high
voltages. Besides, nonlinear electromechanical effects may also occur in
materials with defects, dipole clusters, domains, etc., under moderate fields,
including in ferroelectrics[8, 9, 10]. By definition, the electrostrictive
strain ($x_{ij}$) is proportional to the square of the polarization
($P_{k}P_{l}$) through the electrostrictive fourth-rank tensor $Q_{ijkl}$:

$x_{ij}=Q_{ijkl}\,P_{k}P_{l}$ (1)

We define “linear dielectrics” as insulating materials where the electric
susceptibility ($\chi_{e}$) relating the polarization $P$ to the electric
field $E$ through $P=\epsilon_{0}\chi_{e}\,E$ is constant
($\partial\chi_{e}/\partial E=0$). In these linear dielectrics, a similar
relation (Eq.(2)) relates the strain to the square of the electric field
($E_{k}E_{l}$). The corresponding electrostrictive tensor is then $M_{ijkl}$:

$x_{ij}=M_{ijkl}\,E_{k}E_{l}$ (2)

However, most electrostrictive materials appear to exhibit field-dependent
values of their electrostrictive tensor components rather than constant
values[11, 12, 13, 7] as expected from the definition. The measured
electrostrictive strain is then no longer proportional to the square of the
electric field or of the polarization. Most often, the electrostrictive strain
saturates with increasing electric fields (i.e., the apparent $M$ value
decreases) even though, in some of these cases, the strain nevertheless
remains proportional to the square of the polarization[11] (i.e., the $Q$
values remain constant).

Such a discrepancy between the behavior under the electric field and
polarization can be rationalized as follows. The development of a polarization
inside the dielectric material generates the electrostrictive strain. If the
polarization saturates under large electric fields, then the electric
susceptibility is not linear constant anymore. As a consequence, the increase
in polarization is not anymore proportional to the field increase and the
corresponding strain, while remaining proportional to the saturating
polarization squared, is not proportional anymore to the increasing field
squared beyond the saturation threshold.

The more intriguing case is where the electrostrictive strain saturates as a
function of both the electric field squared and the polarization squared. Such
a case is fundamentally different from the previous one as it also occurs in
linear dielectrics[14], which are, by definition, devoid of polarization
saturation.

Discarding an electric-field dependence on the elastic properties of the
materials, the question of the origin of such saturation remains and is the
topic of this article. We will show that such a saturation is not related to a
field dependence of the electrostrictive coefficients but rather to the
influence of higher-order electromechanical phenomena.

##  2 Nonlinear characterization

Hereafter, we present a method to reliably calculate the value of
electrostrictive tensor components (called “coefficients”) from time-dependent
induced strain measurements. This method enables the fitting of
electrostrictive tensor coefficients that are independent of the electric
field, as is expected from the definition of electrostriction. In addition,
the evolution of the electromechanical strain with the electric field is
explained.

###  2.1 Theoretical analysis

Similar to the definition of nonlinear susceptibilities, nonlinear
electromechanical coupling between strain ($x$) and the electric field ($E$)
or the polarization ($P$) can be expanded using a Taylor expansion[6]:

$x=dE+ME^{2}+M^{(3)}E^{3}+M^{(4)}E^{4}+\ldots$ (3)
$x=gP+QP^{2}+Q^{(3)}P^{3}+Q^{(4)}P^{4}+\ldots$ (4)

where $d$ and $g$ are piezoelectric tensors, $M$ and $Q$ electrostrictive
tensors, the quadratic order of electromechanical interactions. The quantities
$M^{(n)}$ and $Q^{(n)}$ are defined as the $n$-th ($n\geq 3$) order nonlinear
electromechanical tensors. Due to the vector nature of the electric field and
polarization, $M$ and $Q$ are fourth-rank tensors, $M^{(3)}$ and $Q^{(3)}$ are
fifth-rank tensors, and so on.

Centrosymmetric dielectrics exhibit electrostriction but no piezoelectricity
($d$ and $g$ are zero). However, noncentrosymmetric dielectrics can exhibit
both electrostriction and piezoelectricity, the former having been proposed as
an alternative approach to improve the latter[15]. Electrostriction is,
therefore, a second order electromechanical coupling instead of a secondary
effect of piezoelectricity. The contributions to the strain of the higher-
order terms are not necessarily small. They depend on the magnitude of the
fields and on the materials. For example, Park et al. reported that Gd-doped
ceria thin films possess higher-order harmonics[7], which is expressed in
their case through a field-dependent $M$ coefficient. We shall show that
actually, such field-dependence is entirely due to the contribution of higher-
order terms and the $M$ coefficient is constant.

A common electromechanical experiment is to measure the strain of a sample
driven by a mono-frequency electric field:

$E(t)=E_{0}\sin(\omega t)$ (5)

where $E_{0}$ the magnitude of the driving field, $t$ the time, and $\omega$
the angular frequency. Combining Eq.(3) and Eq.(5), the expression of time-
dependent strain induced by a sinusoidal electric field is:

$\begin{split}x(t)=&dE_{0}\sin(\omega t)+ME_{0}^{2}\sin^{2}(\omega t)\\\
&+M^{(3)}E_{0}^{3}\sin^{3}(\omega t)+M^{(4)}E_{0}^{4}\sin^{4}(\omega
t)+\ldots\end{split}$ (6)

In order to fit this time-dependent strain and extract the electrostrictive
coefficients, a model selection procedure could be used to determine the
maximum power to consider, from a trade-off between model complexity and
performance. We show below that this procedure can be advantageously replaced
by using an approach based on the Fourier transform that provides additional
information on the electromechanical response.

The discrete Fourier series ($\\{F(k\omega\\}$) of the time-dependent strain
is defined as:

$F(k\omega)=\sum_{m=0}^{N-1}x(t_{m})\exp{(-i\cdot k\omega\cdot t_{m})}\quad
k=0,\ldots,N-1$ (7)

where $\omega=2\pi/T$ is the excitation angular frequency and $T$ its period,
$N.\delta t=T$ the total acquisition time, $N$ the number of recorded points,
$\delta t$ the sampling interval ($\delta t=t_{m+1}-t_{m}$), and
$x(t_{m})=x(m\delta t)$ the strain recorded at the time $t_{m}$.

The corresponding inverse discrete Fourier series provides the expression of
the strain as a function of the $k$-th harmonics (_i.e._ at $k\omega$) of the
excitation fields:

$x(t)=\frac{1}{N}\sum_{k=0}^{N-1}F(k\omega)\exp{(i\cdot k\omega\cdot t)}$ (8)

The values of the $M$ are sometimes extracted from such a procedure, based on
the Fourier transform of the strain, by only considering the component at
2$\omega$ (_i.e._ for $k$=2 only in Eq.(8)). This can be achieved through the
use of a lock-in amplifier or numerically. We shall show below that the as-
obtained $M$ values (i.e., considering that the strain at $2\omega$ is only
related to the $M$ coefficient) can deviate significantly from the actual
value and exhibit an artificial field dependence, even in purely linear
dielectrics. Indeed, the electrostrictive coefficient $M$ is not the only one
contributing to the strain at $2\omega$.

To show this, we use the power-reduction formulae

$\displaystyle\sin^{2}(\omega t)=$ $\displaystyle\frac{1-\cos(2\omega t)}{2}$
$\displaystyle\sin^{3}(\omega t)=$ $\displaystyle\frac{3\sin(\omega
t)-\sin(3\omega t)}{4}$ $\displaystyle\sin^{4}(\omega t)=$
$\displaystyle\frac{3-4\cos(2\omega t)+\cos(4\omega t)}{8}$ (9)

to express Eq.(6) as a function of multiples of the excitation angular
frequency i.e., in terms of harmonics:

$\displaystyle x(t)=$
$\displaystyle\left[\frac{1}{2}ME_{0}^{2}+\frac{3}{8}M^{(4)}E_{0}^{4}+\ldots\right]+$
(10)
$\displaystyle\left[dE_{0}+\frac{3}{4}M^{(3)}E_{0}^{3}+\ldots\right]\sin(\omega
t)+$
$\displaystyle\left[-\frac{1}{2}ME_{0}^{2}+\frac{1}{2}M^{(4)}E_{0}^{4}+\ldots\right]\cos(2\omega
t)+$
$\displaystyle\left[-\frac{1}{4}M^{(3)}E_{0}^{3}+\ldots\right]\sin(3\omega
t)+$ $\displaystyle\left[\frac{1}{8}M^{(4)}E_{0}^{4}+\ldots\right]\cos(4\omega
t)+\ldots$

where the $\ldots$ stands for higher-order terms.

The third line of Eq.(10) corresponds to the strain occurring at $2\omega$,
that is the second harmonic of the strain. It is a weighted sum of all
electromechanical couplings of even orders ($M\equiv M^{(2)}$, $M^{(4)}$, …)
and does not only consist of electrostriction (_i.e._ $M$ is not the only term
in the prefactor). It is therefore only appropriate to assimilate this pre-
factor to pure electrostriction if the higher-order terms ($M^{(4)}$, …) are
negligible.

To elaborate on the fact that the pre-factor corresponding to the second
harmonic is not only due to electrostriction, Eq.(10) can be re-written as a
sum of cosines (note the sign change of the terms in some pre-factors):

$\displaystyle x(t)=$
$\displaystyle\left[\frac{1}{2}ME_{0}^{2}+\frac{3}{8}M^{(4)}E_{0}^{4}+\ldots\right]+$
(11)
$\displaystyle\left[dE_{0}+\frac{3}{4}M^{(3)}E_{0}^{3}+\ldots\right]\cos\left(\omega
t-\frac{\pi}{2}\right)+$
$\displaystyle\left[\frac{1}{2}ME_{0}^{2}-\frac{1}{2}M^{(4)}E_{0}^{4}-\ldots\right]\cos\left(2\omega
t-\pi\right)+$
$\displaystyle\left[\frac{1}{4}M^{(3)}E_{0}^{3}-\ldots\right]\cos\left(3\omega
t-\frac{3\pi}{2}\right)+$
$\displaystyle\left[\frac{1}{8}M^{(4)}E_{0}^{4}+\ldots\right]\cos\left(4\omega
t-2\pi\right)+\ldots$

The formulation of Eq.(11) is close but not identical to the one of a Fourier
series. To illustrate this, the Fourier series corresponding to Eq.(8) can
also be written as an amplitude-phase form (as $x(t)$ is real, its transform
is hermitian):

$x(t)=\frac{F_{0}}{N}+\frac{2}{N}\sum_{k=1}^{\left\lfloor
N/2\right\rfloor}|F(k\omega)|\,\cos{(k\omega\cdot t-\varphi_{k})}$ (12)

where $\left\lfloor N/2\right\rfloor$ stands for the floor value of $N/2$ to
encompass both odd and even values of $N$, $\varphi_{k}$ is the phase of
$F(k\omega)$ (i.e., $F(k\omega)=|F(k\omega)|e^{i\varphi_{k}}$) and the first
term corresponds to the average value of $x(t)$ over the recorded period.

To further underline the difference between Eq.(12) and Eq.(11), we explicitly
write Eq.(12) up to the fourth harmonic:

$\displaystyle x(t)$ $\displaystyle=$ $\displaystyle\frac{F_{0}}{N}+$ (13)
$\displaystyle\left[\frac{2|F(\omega)|}{N}\right]\cos{(\omega\cdot
t-\varphi_{1})}+$
$\displaystyle\left[\frac{2|F(2\omega)|}{N}\right]\cos{(2\omega\cdot
t-\varphi_{2})}+$
$\displaystyle\left[\frac{2|F(3\omega)|}{N}\right]\cos{(3\omega\cdot
t-\varphi_{3})}+$
$\displaystyle\left[\frac{2|F(4\omega)|}{N}\right]\cos{(4\omega\cdot
t-\varphi_{4})}+\cdots$

Therefore the presence of the phase of the Fourier weights (the
$\varphi_{n}$’s in Eq.(13)) instead of successive multiples of $\pi/2$ (in
Eq.(11)) prevents the identification of the normalized amplitude spectrum of
the Fourier series (the $\left\\{\frac{2|F(k\omega)|}{N}\right\\}$ in Eq.(13))
with the pre-factors (the terms between square brackets) of Eq.(11).

###  2.2 Methodology to calculate electromechanical coefficients

Even though it is not possible to directly calculate the electrostrictive
coefficients (in Eq.11) from the Fourier spectrum of the strain (Eq.13), a
Fourier analysis is nevertheless useful to obtain the electromechanical
coefficients, including the $M$ electrostrictive coefficient. We present below
a reliable method to extract the electromechanical coefficients based on the
measurements of $E(t)$ and $x(t)$ that (1) reproduces the observed saturation
of the electrostrictive strain even in linear dielectrics and that (2)
provides values of the electromechanical coefficients that are field-
independent.

The amplitude spectrum $\left\\{\frac{2|F(k\omega)|}{N}\right\\}$ of the
Fourier transform of $x(t)$ as defined by Eq.(12) is computed first. It serves
as a basis for the choice of the harmonics to consider. The highest harmonic
($h_{m}\omega$) to consider should be set to prevent under- or over-fitting
the data.

It is indeed necessary to only keep harmonics with weights exceeding the
uncertainties of the fast Fourier transform. The substantial aliasing and
spectral leakage errors that discrete Fourier transform may introduce are very
sensitive to measurement noise contaminations[16, 17]. To address this issue,
a method based on the Volterra series was reported for the accurate
measurement of higher-order frequency response functions to identify
structural nonlinearities[18]. Considering the calculation complexity,
measuring entire multi-dimensional higher-order frequency response functions
may not be necessary. For practical purposes, we shall provide a simpler
description to calculate electromechanical coupling coefficients in the
presence of higher-order terms: a histogram of the Fourier weights with a
hundred bin is used to determine the relevant harmonics (see Suppl.
Information). In our measurements, the noise level of the FFT corresponds to 1
or 2% of the normalized Fourier weights. All harmonics corresponding to a
Fourier weight above that threshold are selected.

Once the harmonics are selected, a fit of $x(t)$ is carried out111Eq.14 gather
the electromechanical coupling constants and the powers of the field in order
to have coefficients of similar amplitude, to ease the fitting procedure. The
$x(t)$ curve is fitted rather than the $x(E)$ curve for this reason as well,
in addition to the fact that the losses complicate the fit. with a function
$f(t)$ defined as:

$f(t)=\sum_{\\{h_{n}\\}}f_{n}\sin^{n}(\omega t+\phi_{n})$ (14)

The result of the fit is the set of $\\{f_{n}\\}$ that are the $n$-th
electromechanical coupling coefficients times the corresponding electric field
power: $f_{1}=dE_{0}$, $f_{2}=ME_{0}^{2}$, and $f_{n}=M^{(n)}E_{0}^{n}$
$\forall n\geq 3$. Note that these are not the pre-factors of Eq.(10) but of
Eq.(6). They therefore provide directly (once divided by the $n$-th power of
the electric field) the $n$-th electomechanical coefficient.

##  3 Experimental examples

We illustrate this methodology on classical and non-classical electrostrictors
to demonstrate its ability to both prove field-independent values of the
electrostrictive coefficients and explain the saturation of the
electromechanical strain shown in Figure 1. All measurements are carried out
at 10 Hz and 20∘C on ceramic samples of similar dimensions. The measured
longitudinal displacements induced by the electric field along the same
direction are of several hundred nanometers, far greater than the minimal
measurable value (4 nm) and the measured strains have been averaged over 10
cycles.

###  3.1 Harmonic analysis of classical electrostrictors

Figure 1: Polarization as a function of the electric field ($P$-$E$ loop) of
0.9Pb(Mg_1/3Nb_2/3)O3-0.1PbTiO3 (PMN-10PT) ceramics (a) and substituted
La2Mo2O9 ceramics (d). Electrostrictive responses of PMN-10PT as a function of
the square of electric field amplitude (b) and of the square of polarization
(c). All the measurements are at 10 Hz and 20∘C. The $x$-$E^{2}$ curve of
PMN-10PT (b) saturates while the $x$-$P^{2}$ curve (c) remains linear. In
contrast, substituted La2Mo2O9 are all linear dielectrics (d). As a
consequence, the $x$-$E^{2}$ curve of substituted La2Mo2O9 (e) and $x$-$P^{2}$
(f) curve exhibit a similar tendency: Gd- and (Gd,W)-substituted La2Mo2O9
saturate at the large fields whereas Nd-substituted La2Mo2O9 shows the
opposite behavior.

We call “classical electrostrictors” materials exhibiting electrostrictive
properties described by the empirical proxy[6]. Some of the most common ones
are based on the relaxor ferroelectric lead magnesium niobate (PMN), with a
maximum strain around $0.2\%$[19, 20, 21, 22] for ceramics. Ceramic samples of
0.9Pb(Mg_1/3Nb_2/3)O3-0.1PbTiO3 (PMN-10PT) have therefore been investigated.
The electric-field-induced evolution of polar nano-regions (PNRs) in relaxor
ferroelectrics plays an important role in their electromechanical
properties[23, 24, 25]. In this PMN-10PT, the saturation of the polarization
(Figure 1(a)) causes the saturation of the electrostrictive strain under the
electric field. This is illustrated in Figure 1(b) where the maximum strain
saturates when plotted against the electric field squared, whereas it remains
linear when plotted against the polarization squared (Figure 1(c)). Such
behavior leads to an apparent non-constant, field-dependent $M$ coefficient,
whereas the $Q$ coefficient remains constant as expected.

The curve $x$-$E^{2}$ (Figure 1(b)) for PMN-10PT ceramics may be thought of as
consisting of two parts: a linear part below about 1 MV/m and a non-linear
part under higher electric field. Li et al. explained similar results in PMN
crystals by separating the response into such low- and high-field regions[26].
In the high-field region, the saturation of the electromechanical strain was
explained by the fact that all the PNRs have been switched to a stable state,
causing the saturation of both polarization and electrostriction. This
explanation is indeed consistent with Figure 1(c), showing that the curve
$x$-$P^{2}$ is approximately linear rather than saturated. In relaxor
ferroelectrics as well as in nonlinear dielectrics, the saturation of the
polarization under the electric field is indeed usually the primary reason for
the saturation of the electromechanical strain. This is the reason why
electrostriction in ferroelectrics is preferably assessed by $Q$ coefficients.
Sundar et al. demonstrated in addition that the $Q$ coefficient is fairly
independent of temperatures in PMN ceramics[27]. We shall however show below
that this explanation omits the contribution of higher-order couplings and
that such separation in low- and high-field parts is actually not necessary.

Figure 2: Nonlinear strain responses of PMN-10PT ceramic as a function of
electric field (a) and time (b) in the low-field region (0.5 MV/m) and high-
field region (2 MV/m) (d and e, respectively). Amplitude spectra from the
Fourier transform represent the weight of the $n$-th harmonics for low-(c) and
high-fields(f). The dashed lines in (b) and (e) are the fitted strains using
Eq.14. The hollow symbols in (c) and (f) represent the discarded harmonics,
and the fit is carried out considering only the contribution from the
harmonics represented by solid symbols. (g) The “apparent” $M$ coefficient,
calculated from the deformation at 2$\omega$ (therefore including higher-order
contributions) decreases as a function of the electric field whereas (h) the
$M$ coefficient (corresponding to pure electrostriction) remains constant at
about $1.76(\pm 0.15)\times 10^{-15}\,\rm m^{2}V^{-2}$ for all field
amplitudes.

In order to test Li’s description based on low-field and high-field regions,
ac electric fields between 0.5 MV/m a value approximatively equal to the
coercive field (Figures 2(a-c)) and 2 MV/m (Figures 2(d-f)) are applied to the
PMN-10PT sample of Figures 1(a-c). Figures 2(a) and (d) show the strain that
develops under the minimum and maximum amplitudes of electric fields
considered, indicating that electrostriction is the predominant
electromechanical coupling in PMN-10PT ceramics regardless of the field
amplitude, as is expected. Correspondingly, Figures 2(b) and (e) show that the
strain is positive and evolves, as a first approximation, at twice the
frequency of the electric field (the strain is measured over one period of the
electric field). The normalized amplitude spectra calculated by Fourier
analysis for each field amplitude are presented in Figures 2(c) and (f). At
low field, Figure 2(c) indicates that the second harmonic predominates (as
expected) but also that the first (piezoelectric) and fourth harmonic
($M^{(4)}$) contributions are not negligible. The dashed line in Figure 2(b)
represents the result of a sinusoidal fit of $x(t)$ using Eq.(14) with linear,
quadratic, and quartic terms ($n=1,2,4$) and is a satisfactory reproduction of
the measured strain with little residuals.

In the high-field region, the comparison of Figure 2(c) and (f) shows that the
normalized intensities of higher-order harmonics increase. This is an expected
consequence of increasing the maximum field in a non-linear system. Only
keeping the second, fourth, sixth, eighth, and tenth harmonics in Figure 2(f),
this saturation curve can also be satisfactorily reproduced by fitting $x(t)$
Eq.(14) with $n$=2,4,6,8, and 10 (see Supplementary Information, Figure S1).

The above results have several consequences: (1) the calculated $M$ in the
low-field region is $1.58(\pm 0.04)\times 10^{-15}\,\rm m^{2}V^{-2}$,
consistent with the value calculated in the high-field region: $1.93(\pm
0.03)\times 10^{-15}\,\rm m^{2}V^{-2}$, when the higher-order terms are taken
into account. Hence, the $M$ remains constant (within less than 20%)
regardless of the amplitude of the driving field, as shown in Figure 2(h).
This is the expected behavior of a “constant” (see Supplementary Information
about the evolution of $M^{(n)}$, Figure S2). (2) The low-field and high-field
classification does not even ensure that PMN-10PT is purely electrostrictive
(i.e., with no higher-order electromechanical terms), as even in the low-field
region, the contribution of the quartic order interaction is not negligible.
If only the strain at 2$\omega$ is considered (by filtering the Fourier
spectrum to consider only this harmonic) then the $M$ coefficient value is
1.09$\times$10-15 m2V-2 at low field and 0.37$\times$10-15 m2V-2 at high
field, as shown in Figure 2(g). Such apparent field dependence of the $M$
coefficient is therefore due to the omission of higher-order coupling
coefficients.

Therefore, in PMN-10PT, the saturation of the polarization causes the
saturation of the maximum electromechanical strain. Such a saturation should
not be attributed to a field dependence of the $M$. In addition, calculating
the $M$ coefficient from a filter centered on the $2\omega$ frequency leads to
erroneous values.

From a microscopic perspective, the active PNRs behave as “seeds” to induce
the macro-domain from the dielectric matrix and then facilitate the switching
of the macro-domains. The precise description of this dynamic process reflects
the nonlinear electrostriction-based nature of relaxor ferroelectrics[28].

###  3.2 Harmonic analysis of non-classical electrostrictors

On the contrary to “classical electrostrictors”, some materials do not follow
the empirical proxy based on elastic and dielectric properties[6]. When the
deviation exceeds one order of magnitude, these materials have been defined as
“giant” or “non-classical” electrostrictors[4]. Yavo et al. have investigated
the saturation of electrostriction in doped ceria ceramics, the first non-
classical electrostrictor reported, pointing to charge trapping, voltage
redistribution at the grain boundaries, and electrodes as origins of the
observed saturation[29, 13, 12, 5]. La2Mo2O9-based materials (LAMOX) have also
been reported to exhibit non-classical electrostriction [30, 31]. Hereafter,
we report results on ceramics of three of these LAMOX materials (structural
and microstructural information can be found in the Supplementary Information,
Figure S3 and S4).

In contrast to the classical electrostrictor case, the maximum strain of
substituted La2Mo2O9 ceramics, i.e. for non-classical electrostrictors,
exhibits a similar non-linear evolution when plotted as a function of $E^{2}$
and $P^{2}$ (see Figure 1(e) and (f)). This is consistent with the linear
dielectric character (Figure 1(d)) of LAMOX ceramics that allows to link the
two electrostrictive coefficients to be made through $M=Q\varepsilon^{2}$ with
$\varepsilon$ the permittivity.

La2Mo2O9 and W-substituted La2Mo2O9 ceramics were reported to exhibit non-
classical electrostrictive effects that can be enhanced at elevated
temperatures[30]. It was further revealed that oxygen vacancy configurations
and hopping dynamics cause the distortion of the local structure without
measurable bulk polarization[30]. Here, we shall extend the study of this
family of materials by considering other substitutions.

Note that the curves of $x$-$E^{2}$ (Figure 1(e)) and $x$-$P^{2}$ (Figure
1(f)) of Gd-substituted (La_1.8Mo2Gd_0.2O9) and (Gd,W)-substituted La2Mo2O9
(La_1.7Mo_0.8Gd_0.3W_1.2O9) saturate, whereas the Nd-substituted La2Mo2O9
(La_1.4Mo2Nd_0.6O9) curves shows the opposite tendency i.e., they tend to
_increase_ under increasing electric fields. Clearly, the explanation posited
for classical electrostrictors (that the saturation of the polarization causes
the saturation of the strain under electric field) can not be applied to non-
classical (“giant”) electrostrictors such as LAMOX as they remain linear
dielectrics (see Figure 1(d)). Even though these “giant” electrostrictors have
been studied from different angles[32, 33, 34], only hypotheses have been put
forward to explain their strain saturation. We show below that such behavior
is explained when higher-order electromechanical coefficients are considered,
thereby providing a unified description of the electromechanical strain
saturation in both classical and non-classical electrostrictors.

Figure 3: Negative electrostrictive response of Gd-substituted La2Mo2O9
ceramic as a function of the electric field (a) and time (b). The sinusoidal
fit (dashed line in (b)) using Eq.14 reproduces the measured strain (blue
line) by only considering the orders corresponding to harmonics which are
identified from the normalized intensity of the fast Fourier transform of
$x(t)$ (c). The opposite signs of $M=-1.89\times 10^{-17}\,\rm m^{2}V^{-2}$
and $M^{(4)}$=$+3.27\times 10^{-31}\,\rm m^{4}V^{-4}$ induce the saturation of
the electrostrictive strain.

Figure 3(a) shows the field-induced strain in Gd-substituted La2Mo2O9
indicative of a predominant negative electrostriction. Figure 3(b) presents
the time-dependence of the strain, measured over one field period and
exhibiting at first glance an evolution at twice the field frequency. Both of
these curves appear to point to a predominant electrostrictive behavior.
However, the amplitude spectrum (Figure 3(c)) shows that, in addition to the
(expected) predominant second harmonic, the first and fourth harmonics are not
negligible. The fitting of the $x(t)$ curve (represented by the dash line in
Figure 3(b)) provides a negative value for $M=-1.89(\pm 0.02)\times
10^{-17}\,\rm m^{2}V^{-2}$ and a positive value for $M^{(4)}$=$+3.27(\pm
0.14)\times 10^{-31}\,\rm m^{4}V^{-4}$. Therefore, as the field increases, the
positive $M^{(4)}E_{0}^{4}$ term will compensate more and more for the
negative $ME_{0}^{2}$ term. The saturation of the maximum strain in Gd-
substituted La2Mo2O9 is therefore due to the opposite signs of $M$ and
$M^{(4)}$.

Another consequence of taking into account the fourth harmonic is that it
leads to a value that differs from the one that would be obtained if only the
response at $2\omega$ was considered ($-0.96\times 10^{-17}\,\rm
m^{2}V^{-2}$). This underlines, as was the case with classical
electrostrictors, the importance of taking into account higher-order
electromechanical terms to accurately evaluate the values of electrostrictive
coefficients. It is worth pointing out that the non-negligible first harmonic
does not necessarily correspond to piezoelectricity. Indeed, the last measured
strain point differs from the first measured point (see Figure 3(a), which can
also be seen in the residuals in Figure 3(b)). This can be due to various
phenomena, e.g., to the damping behavior of defect dipoles or to the initial
alignment of dipoles during the first application of the electric field. This
leads to a non-zero overall slope of the strain as a function of time that
contributes to the non-zero first harmonic.

Figure 4: Negative electrostrictive response of (Gd-W)-substituted La2Mo2O9
ceramic as a function of the electric field (a) and time (b). The sinusoidal
fit (dashed line in (b)) using Eq.14 reproduces the measured strain (blue
line) by only considering the orders corresponding to harmonics which are
identified from the normalized intensity of the fast Fourier transform of
$x(t)$ (c). The opposite signs of $M=-0.35\times 10^{-17}\,\rm m^{2}V^{-2}$
and $M^{(4)}=+5.53\times 10^{-32}\,\rm m^{4}V^{-4}$ induce the saturation of
the electrostrictive strain.

If the additional addition of tungsten in Gd-substituted La2Mo2O9 enables to
reduce the electrostrictive losses (the strain vs electric field curve is less
hysteretic than without tungsten, compare Figure 3(a) and Figure 4(a)), it
comes at the cost of decreasing significantly (5-fold reduction) the maximum
strain. The $M$ coefficient of (Gd,W)-substituted La2Mo2O9 ($M=-0.35(\pm
0.01)\times 10^{-17}\,\rm m^{2}V^{-2}$), considering the first, second, and
fourth harmonic is therefore lower than the one of Gd-substituted La2Mo2O9. As
in the case of Gd-substituted La2Mo2O9, the apparent saturation of
electrostriction remains due to the opposite signs of $M$ and $M^{(4)}$
($+5.53(\pm 0.56)\times 10^{-32}\,\rm m^{4}V^{-4}$).

We now turn to Nd-substituted La2Mo2O9. This material exhibits a positive
electrostrictive strain (see Figure 5(a)) with an electrostrictive coefficient
$M=+0.49(\pm 0.01)\times 10^{-17}\,\rm m^{2}V^{-2}$ similar in amplitude but
opposite in sign to the one of (Gd,W)-substituted La2Mo2O9. In addition to the
second harmonic, the first, and third harmonics are considered from the
amplitude spectrum (Figure 5(c)) in the fitting of $x(t)$ (Figure 5(b)) using
Eq.14 when $n=1,2,3$. In contrast to the previous cases, the quartic order
term is negligible, while the $M^{(3)}$ coefficient ($+5.36(\pm 1.27)\times
10^{-26}\,\rm m^{3}V^{-3}$) exhibits a positive sign. As a consequence, rather
than compensating each other and causing the electromechanical strain to
saturate as the field increases, the electromechanical strain of Nd-
substituted La2Mo2O9 _increases_ as the field amplitude increases. The modest
relative weight of the harmonics considered compared to the amplitude of the
second harmonic explains the modest positive deviation from linearity over the
investigated field range. In addition, as $M^{(3)}$ does not contribute to the
signal at 2$\omega$ (cf. Eq.10), the value of $M$ one would get considering
only the second harmonic (0.44$\times$10-17m2V-2) is close (within 10%) to the
actual one. However, not considering this higher-order coupling prevents from
explaining the increase of the electromechanical strain when the field
increases.

Figure 5: Positive electrostrictive response of Nd-substituted La2Mo2O9
ceramic as a function of the electric field (a) and time (b). The sinusoidal
fit (dashed line in (b)) using Eq.14 reproduces the measured strain (blue
line) by only considering the orders corresponding to harmonics which are
identified from the normalized intensity of the fast Fourier transform of
$x(t)$ (c). The similar signs of $M=+0.49\times 10^{-17}\,\rm m^{2}V^{-2}$ and
$M^{(3)}=+5.36(\pm 1.27)\times 10^{-26}\,\rm m^{3}V^{-3}$ induce the increase
of the electrostrictive strain when the field amplitude increases.

The conclusion to draw from these substituted La2Mo2O9 is that substitution
elements can tune the amplitude, hysteresis, and saturation behavior of the
electromechanical properties in non-classical electrostrictors, through the
contributions of higher-order electromechanical coupling coefficients.
Considering higher-order harmonics based on the amplitude spectra of their
Fourier transform enables in non-classical electrostrictors as in classical
electrostrictors to calculate a field-independent and more accurate values of
the electrostrictive coefficients. We anticipate that microstructural analyses
of these substituted La2Mo2O9 will enable to correlate the signs of the
higher-order coupling coefficients to microstructural features.

##  4 Conclusion

The study of classical and non-classical electrostrictive materials requires
the consideration of higher-order electromechanical terms, even at low fields,
to calculate the actual values of electrostrictive coefficients. The apparent
field-dependence of electrostrictive coefficients calculated either from a
quadratic fit of the strain versus electric field curves or after having
filtered only the second harmonic may explain the dispersion of
electrostrictive coefficients reported in the literature. Indeed, such a
procedure does not enable the determination of an accurate value of the
electrostrictive coefficients. In addition, the signs of the higher-order
coefficients enable predictions to be made about the behavior of the maximum
strain under increasing electric field. The saturation (or increase) of the
electromechanical strain under electric field is explained by the increasing
contribution of higher-order couplings in the electromechanical response.

We propose a scheme to determine the actual values of the electrostrictive
coefficients and predict the evolution of electromechanical strain under
electric field: first, the amplitude spectrum of the time-dependent strain
curve is calculated from its Fourier transform, then the harmonics to consider
are chosen based on the histogram of this amplitude spectrum, and finally are
used to carry out a fit of the time evolution of the strain with a polynomial
function of $\sin(\omega t)$. This provides values of electrostrictive
coefficients that do not depend on the amplitude of the applied electric
field. In addition, from the sign of the higher-order electromechanical
coefficients, it is possible to determine whether they will compensate each
other (if they have opposite signs) or add to each other (if they have the
same sign) and therefore predict whether the strain will saturate or increase
when the field increases.

Acknowledgements This work was partially financed by ANR-19-ASTR-0024-01 and
ANR-20-CE08-0012-1 grants

## References

  * [1] R. Korobko, A. Patlolla, A. Kossoy, E. Wachtel, H. L. Tuller, A. I. Frenkel, I. Lubomirsky,  _Advanced Materials_ 2012, _24_ , 43 5857.
  * [2] N. Yavo, A. D. Smith, O. Yeheskel, S. Cohen, R. Korobko, E. Wachtel, P. R. Slater, I. Lubomirsky,  _Advanced Functional Materials_ 2016, _26_ , 7 1138.
  * [3] E. Makagon, O. Kraynis, R. Merkle, J. Maier, I. Lubomirsky,  _Advanced functional materials_ 2021, _31_ , 50 2104188.
  * [4] J. Yu, P.-E. Janolin,  _Journal of Applied Physics_ 2022, _131_ , 17.
  * [5] M. Varenik, B. Xu, J. Li, E. Gaver, E. Wachtel, D. Ehre, P. K. Routh, S. Khodorov, A. I. Frenkel, Y. Qi, et al.,  _Nature Communications_ 2023, _14_ , 1 7371.
  * [6] R. Newnham, V. Sundar, R. Yimnirun, J. Su, Q. Zhang,  _The Journal of Physical Chemistry B_ 1997, _101_ , 48 10141.
  * [7] D.-S. Park, M. Hadad, L. Riemer, R. Ignatans, D. Spirito, V. Esposito, V. Tileli, N. Gauquelin, D. Chezganov, D. Jannis, et al.,  _Science_ 2022, _375_ , 6581 653.
  * [8] L. M. Riemer, L. Jin, H. Uršič, M. Otonicar, T. Rojac, D. Damjanovic,  _Journal of Applied Physics_ 2021, _129_ , 5.
  * [9] A. Grigoriev, R. Sichel, H. N. Lee, E. C. Landahl, B. Adams, E. M. Dufresne, P. G. Evans,  _Physical review letters_ 2008, _100_ , 2 027604.
  * [10] D. Damjanovic,  _Reports on progress in physics_ 1998, _61_ , 9 1267.
  * [11] L. Cross, S. Jang, R. Newnham, S. Nomura, K. Uchino,  _Ferroelectrics_ 1980, _23_ , 1 187.
  * [12] M. Varenik, J. C. Nino, E. Wachtel, S. Kim, O. Yeheskel, N. Yavo, I. Lubomirsky,  _ACS applied materials & interfaces_ 2020, _12_ , 35 39381.
  * [13] A. Kabir, S. Santucci, N. Van Nong, M. Varenik, I. Lubomirsky, R. Nigon, P. Muralt, V. Esposito,  _Acta Materialia_ 2019, _174_ 53\.
  * [14] K. Uchino,  _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ 2022, _69_ , 11 3013.
  * [15] F. Li, L. Jin, Z. Xu, S. Zhang,  _Applied Physics Reviews_ 2014, _1_ , 1.
  * [16] Y. Ma,  _IEEE transactions on signal processing_ 1997, _45_ , 6 1641.
  * [17] H. Dishan,  _Mechanical Systems and Signal Processing_ 1995, _9_ , 2 113.
  * [18] R. Lin, T. Ng,  _ISA transactions_ 2018, _81_ 270\.
  * [19] S.-E. Park, T. R. Shrout,  _Journal of applied physics_ 1997, _82_ , 4 1804.
  * [20] S.-F. Liu, S.-E. Park, T. R. Shrout, L. Eric Cross,  _Journal of applied physics_ 1999, _85_ , 5 2810.
  * [21] A. Bokov, Z.-G. Ye,  _Journal of materials science_ 2006, _41_ 31\.
  * [22] N. Donnelly, G. Catalan, C. Morros, R. Bowman, J. Gregg,  _Journal of applied physics_ 2003, _93_ , 12 9924.
  * [23] G. Burns, F. Dacol,  _Physical Review B_ 1983, _28_ , 5 2527.
  * [24] G. Shirane, G. Xu, P. Gehring,  _Ferroelectrics_ 2005, _321_ , 1 7.
  * [25] X. Shi, J. Wang, J. Xu, X. Cheng, H. Huang,  _Acta Materialia_ 2022, _237_ 118147\.
  * [26] F. Li, Z. Xu, S. Zhang,  _Applied Physics Letters_ 2014, _105_ , 12.
  * [27] V. Sundar, R. Newnham,  _Ferroelectrics_ 1992, _135_ , 1 431.
  * [28] R.-X. Wang, J. Zhang, L.-M. Zheng, L. Sun, Z.-L. Luo, S.-T. Zhang, B. Yang,  _Ceramics International_ 2018, _44_ , 15 18084.
  * [29] N. Yavo, O. Yeheskel, E. Wachtel, D. Ehre, A. I. Frenkel, I. Lubomirsky,  _Acta Materialia_ 2018, _144_ 411\.
  * [30] Q. Li, T. Lu, J. Schiemer, N. Laanait, N. Balke, Z. Zhang, Y. Ren, M. A. Carpenter, H. Wen, J. Li, et al.,  _Physical Review Materials_ 2018, _2_ , 4 041403.
  * [31] J. Yu,  Ph.D. thesis, Université Paris Saclay, 2022.
  * [32] T. Das, J. D. Nicholas, B. W. Sheldon, Y. Qi,  _Physical Chemistry Chemical Physics_ 2018, _20_ , 22 15293.
  * [33] M. Asher, O. Diéguez,  _APL Materials_ 2019, _7_ , 4.
  * [34] D. P. Trujillo, A. Gurung, J. Yu, S. K. Nayak, S. P. Alpay, P.-E. Janolin,  _npj Computational Materials_ 2022, _8_ , 1 251.


Further author information: (Send correspondence to Jacques Kluska)
Jacques Kluska: E-mail<EMAIL_ADDRESS>Telephone: (+32) 16 32 84
51

# VLTI images of circumbinary disks around evolved stars

Jacques Kluska  KULeuven, Celestijnenlaan 200D, 3001 Leuven (Belgium)  Rik
Claes  KULeuven, Celestijnenlaan 200D, 3001 Leuven (Belgium)  Akke Corporaal
KULeuven, Celestijnenlaan 200D, 3001 Leuven (Belgium)  Hans Van Winckel
KULeuven, Celestijnenlaan 200D, 3001 Leuven (Belgium)  Javier Alcolea
Observatorio Astronómico Nacional (Spain)  Narsireddy Anugu  University of
Arizone (United States)  Jean-Philippe Berger  Institut de Planétologie et
d’Astrophysique de Grenoble (France)  Dylan Bollen  KULeuven, Celestijnenlaan
200D, 3001 Leuven (Belgium)  Valentin Bujarrabal  Observatorio Astronómico
Nacional (Spain)  Robert Izzard  University of Surrey (UK)  Devika Kamath
Macquarie University (Australia)  Stefan Kraus  University of Exeter (UK)
Jean-Baptiste Le Bouquin  Institut de Planétologie et d’Astrophysique de
Grenoble (France)  Michiel Min  SRON Netherlands Institute for Space Research
(Netherlands)  John D. Monnier  Univ. of Michigan (United States)  Hans
Olofsson  Chalmers University of Technology (Sweden)

###### Abstract

The new generation of VLTI instruments (GRAVITY, MATISSE) aims to produce
routinely interferometric images to uncover the morphological complexity of
different objects at high angular resolution. Image reconstruction is,
however, not a fully automated process. Here we focus on a specific science
case, namely the complex circumbinary environments of a subset of evolved
binaries, for which interferometric imaging provides the spatial resolution
required to resolve the immediate circumbinary environment. Indeed, many
binaries where the main star is in the post-asymptotic giant branch (post-AGB)
phase are surrounded by circumbinary disks. Those disks were first inferred
from the infrared excess produced by dust. Snapshot interferometric
observations in the infrared confirmed disk-like morphology and revealed high
spatial complexity of the emission that the use of geometrical models could
not recover without being strongly biased. Arguably, the most convincing proof
of the disk-like shape of the circumbinary environment came from the first
interferometric image of such a system (IRAS08544-4431) using the PIONIER
instrument at the VLTI. This image was obtained using the SPARCO image
reconstruction approach that enables to subtract a model of a component of the
image and reconstruct an image of its environment only. In the case of
IRAS08544-4431, the model involved a binary and the image of the remaining
signal revealed several unexpected features. Then, a second image revealed a
different but also complex circumstellar morphology around HD101584 that was
well studied by ALMA. To exploit the VLTI imaging capability to understand
these targets, we started a large programme at the VLTI to image post-AGB
binary systems using both PIONIER and GRAVITY instruments.

###### keywords:

Image reconstruction, infrared interferometry, binaries, VLTI, PIONIER

##  1 INTRODUCTION

A significant fraction of stars are binaries and binary research constitutes a
main domain of stellar astrophysics. Binary evolution gives birth to diverse
phenomena such as thermonuclear novae, supernovae type Ia, sub-luminous
supernovae, gravitational waves and objects such as subdwarf B-stars, barium
stars, cataclysmic variables, and asymmetric planetary nebulae. However,
binary evolution is not yet well understood. Here, we focus on a peculiar
stage of binary evolution that are the post-asymptotic giant branch (post-AGB)
binaries. During the preceding phase, the AGB phase, the star has strong winds
that expel its atmosphere. In the case of a close binary (with a separation
less than few astronomical units) this interaction produces a circumbinary
disk that is observable in the post-AGB phase through the infrared excess in
the spectral energy distribution[1, 2] reminiscent of young stellar objects.
As the binary nature is difficult to reveal during the AGB phase, because of
stellar pulsations and atmosphere motions that blur the radial velocity signal
but also because of the dense dusty surroundings of the star, the post-AGB
phase is ideal to probe binarity and the system surroundings as the atmosphere
of the evolved star is almost all expelled and the direct environment is not
obscuring the stars anymore. We can, therefore, detect and characterise the
orbits of post-AGB binary systems[3]. However, their periods and
eccentricities do not correspond to the expected ones from binary interaction
models[4].

It is now well known that binarity and the presence of a circumbinary disk are
tightly linked[5]. Disk-binary interaction can therefore be the key to solve
the discrepancy between theory and observations. Those disks were further
studied by millimetre interferometry showing their stability as they are in
Keplerian rotation[6, 7, 8, 9, 10]. However, the angular scale achieved with
those observations do not allow to probe the inner disk regions where the
disk-binary interaction processes happen.

To explore this, one needs to study the inner disk structure. As the distances
to those objects are relatively large ($\sim$1 kpc), those disks appear
compact on the sky and the use infrared interferometry is necessary to achieve
the needed angular resolution power to spatially resolve them. In this paper
we show how infrared interferometry offers unique capabilities to constrain
the properties of the systems and how the interferometric imaging technique
contributed to advance our knowledge on these peculiar objects.

In this paper, we will describe the first snapshot-like studies and surveys
using infrared interferometers with sparse (u, v)-coverage (Sect. 2). Then we
will show what interferometric imaging brings to the understanding of these
systems following two examples: IRAS08544-4431 (Sect. 3) and HD101584 (Sect.
4). We will finish by presenting INSPIRING: an interferometric imaging large
programme at the VLTI on post-AGB binaries (Sect. 5).

##  2 Interferometric snapshots

First infrared interferometric snapshot observations were performed using the
MIDI and AMBER instruments on the VLTI strengthening the idea that the
circumbinary environment has a compact disk-like structure[11, 12]. Further
observations using VLTI and CHARA interferometric instruments confirmed it as
well and allowed a successful in-depth modelling of both the SED and the
interferometric data using radiative transfer models of protoplanetary
disks[13, 14, 15]. One of those models based on AMBER and MIDI data allowed to
localise the PAH emission observed in a unique target, HR4049, to be in the
outflow[16].

In order to confirm the first conclusions from observations of individual
targets, interferometric surveys were carried out subsequently. First, in the
mid-infrared[17] using MIDI, a survey of 19 objects have showed a strong
similarity between the observations and a synthetic population of models of
protoplanetary disks, demonstrating that the disk morphology is best suited to
reproduce both the interferometric observations and the SED for a larger
sample of targets. Moreover, a similar mid-infrared size-color relation was
found for those disks as for disks around young stellar objects[18].

A similar snapshot survey was obtained in the near-infrared on a sample of 23
sources[19] using the PIONIER[20] instrument on the VLTI. Using a set of
geometrical models with increasing complexity it was possible to show that the
near-infrared size/stellar luminosity relation found for young stellar
objects[21, 22], can be extrapolated to the high-luminosities that are typical
for post-AGB stars ($10^{3}-10^{4}$ L⊙). This demonstrates that the inner
dusty rims of disks around post-AGB binaries are not ruled by dynamical
truncation due to the inner binary but by dust sublimation. Those near-
infrared sizes of the environment are, however, slightly larger than in
protoplanetary disks likely because of different dust chemistry and/or gas
density in the dust sublimation region. This is corroborated by the slightly
lower environmental temperatures measured using the spectral channels of the
PIONIER instrument. The other main conclusion of this survey was that complex
geometric models (more than 13 parameters) are favoured for most ($>$55%)
targets. This suggests a high intrinsic complexity of the target morphologies
in the near-infrared. There is, therefore, a risk of being model biased when
trying to reproduce the complex interferometric signal provided by 4-telescope
infrared interferometric instruments. With the advent of instruments being
able to combine more and more telescopes in the infrared, interferometric
image reconstruction techniques are perfect to solve this problem and to
reveal the origin of the complexity of these targets.

##  3 First image: IRAS08544-4431

Figure 1: Image reconstructions of IRAS08544-4431 from VLTI/PIONIER data[23,
24]. Left: monochromatic image reconstruction. The flux level was cut to
better see the environment of the star. Center: SPARCO/MiRA image
reconstruction subtracting a single star with diameter 0.5 mas,
$f^{*}_{0}$=61.0% and $d_{\mathrm{env}}$=0.182. This image has revealed
emission coming from the secondary. Right: SPARCO/MiRA image reconstruction
subtracting the primary star and contribution from the secondary, revealing
the circumbinary disk structure. The cyan disk indicates the position and size
of the primary star. The yellow cross indicates the location of the emission
from the secondary star. The beam size is indicated in the bottom-right
corner.

The first interferometric image in the infrared was performed using the
VLTI/PIONIER instrument on IRAS08544-4431 and revealed a very complex
structure[23] as can be seen in Fig. 1. To be able to retrieve the image of
the circumstellar environment, the Semi-Parametric Approach for Reconstruction
of Chromatic Objects (SPARCO) proved to be crucial, as it assigns a
geometrical model and a spectrum to the central star and removes it from the
image, allowing a dramatic increase of the quality of the reconstructed image
of the environment[25]. This technique was successfully applied to image the
environments of young stellar objects[26, 27] and is particularly adapted to
circumstellar environments.

As a first step we have modelled the star as uniform disk with a diameter of
0.5 mas as predicted from the SED modeling. We also fitted the stellar-to-
total flux ratio ($f^{*}_{0}$) to be $\sim$61% of the total flux at 1.65$\mu$m
and the spectral index ($d_{\mathrm{env}}$) to be $\sim$0.182 (defined as
$d_{\mathrm{env}}=\frac{\mathrm{d}\log F_{\lambda}}{\mathrm{d}\log\lambda}$,
and that corresponds in that case to a temperature of $\sim$1700 K). In the
resulting image we can see a clear ring but also a point source emission close
to the primary. This emission is likely associated with the secondary star. It
was unexpected as the secondary is star is likely a main-sequence star and,
hence, the contrast should be such that it is not detectable in our
observations. As we explain it later on, we interpreted this emission as
likely coming from a circum-secondary accretion disk.

We therefore used SPARCO to subtract the emission from the position of the
secondary as well. We used a geometrical model to find the best fitting
positions and flux contribution of the secondary and use them in the SPARCO
image reconstruction. The resulting image provides a better quality
reconstruction of the inner disk rim. The disk is seen almost pole-on and has
complex azimuthal structures showing two maxima and two minima.

We then modelled the data using the radiative transfer code MCMax that was
developped for modeling protoplanetary disks[24]. The model is able to
reproduce both the SED and the squared visibilities from PIONIER. However,
given that the model is axi-symmetric, it was not possible to reproduce the
closure phase signal. This confirms that the structure we see in the image
reconstruction has an astrophysical origin. Moreover, there is a large amount
($\sim$15%) of over-resolved flux (the visibility level at the shortest
baseline is below one) that could not be reproduced by the radiative transfer
model alone.

Finally, this imaging experience allowed us to build an archetype of the inner
regions of these post-AGB binaries with the following building blocks:

  * •

the post-AGB star that emits $\sim$60% of the total flux in $H$-band.

  * •

the inner rim of the circumbinary disk is at the dust sublimation radius and
is azimuthally perturbed

  * •

an unexpected emission from the secondary main-sequence star. This secondary
star should not be detectable given the expected contrast with the post-AGB
star. We hypothesised that the emission comes from the circum-secondary
accretion disk. This is supported by the systematic detection of jets
originating from the secondary star in Hα time-series spectra[28, 29, 30, 31].

  * •

A strong contribution from an over-resolved emission is detected. The origin
of this emission is unknown yet and cannot be reproduced by radiative transfer
model including a disk alone.

Interferometric imaging was crucial to reveal such a detailed picture of this
complex system which is an archetype for a large number of evolved binaries.
However, it is yet unclear if IRAS08544-4431 is a typical example of post-AGB
binary systems. One needs to study other targets of this class.

##  4 Revealing a complex environment: imaging HD101584

Figure 2: Two physical scales probed by ALMA and PIONIER. Left: ALMA image in
CO(2–1) with a resolution of 85 mas. Right: VLTI/PIONIER image with a 2 mas
resolution. The star, that is subtracted from the image using SPARCO, is
indicated by the green cross. The beam of the VLTI observations is in the
bottom-right corner.

Figure 3: VLTI/PIONIER image reconstructions of HD101584 using
MiRA/SPARCO[32]. Left: image with $f^{*}_{0}$=24.4% and $d_{\mathrm{env}}$ =
2.45. Center and right: image with $f^{*}_{0}$=28.0% and $d_{\mathrm{env}}$ =
-1.75. On the right image the contours indicate the 5-$\sigma$ significance of
the emission as estimated by the bootstrapping.

Another infrared interferometric imaging campaign was performed on HD101584
that has been extensively studied with the Atacama Large
Millimeter/submillimeter Array (ALMA)[33, 34, 35]. These observations revealed
a bipolar environment as traced by the CO gas (see Fig. 2). An hourglass
structure is dominating the large scale morphology with high-velocity
outflows. Close to the star, there is a disk-like structure seen both in CO
and in the continuum. An environment with such a symmetry axis is likely to be
the product of a strong binary interaction. However, photometric[36] and
spectroscopic[37] time-series do find a periodic signal but with different and
contradicting periods (218 days and 144 days resp.) and the detections do not
seem to be significant enough. Today’s binary status of the target is
therefore unclear, but it is likely that, if not binary today, the object has
been a binary in the past and merged.

The observations of the inner regions of this target were made with
VLTI/PIONIER. We used the SPARCO technique again to recover the image from the
interferometric dataset. Here, we subtracted a single point-source and
retrieved a stellar-to-total flux ratio of $f^{*}_{0}$ of 24.4% and
$d_{\mathrm{env}}$ of 2.45. The image shows a complex morphology that can be
divided into an inner and an outer morphology (see Fig. 3). The inner
morphology appears as a ring. The outer morphology seems also to form a ring-
like feature even though the emission is not azimuthally continuous. This
second ring seems to be shifted compared to the first ring. No secondary is
detected. It could be due to the lack of angular resolution or because the
target is a merger [32]. The rest of the emission is not significant as
estimated from bootstrapping.

This first image inspired the geometrical models that include two rings, a
central star and a background emission. The geometrical models have shown that
the inner ring has a temperature of 1540$\pm$10 K while the outer one is
colder with a temperature of 1014$\pm$10 K. This geometrical model confirmed
the shift between the two rings with a fitted value of 0.57$\pm$0.01 mas. It
is interesting to note that the shift is in the direction of the minor axis as
it would be the case for an emission above the disk mid-plane. The stellar-to-
total flux ratio obtained from the geometrical model is 28%. As a cross-check,
we have used this stellar-to-total flux ratio to produce a new image
reconstruction. We can see that the regions around the star are depleted from
emission as expected with a higher $f^{*}_{0}$ (see Fig. 3).

Several scenarios could explain the results from interferometric imaging and
geometrical model fitting (see the discussion in [32]) but here is what we
consider the most likely scenario. The inner ring size (4.1 au in radius) is
compatible with the dust sublimation rim of the circumstellar disk. The outer
ring does not lie in the plane of the disk and is above at a height of 1.6 au.
We speculate that the outer ring could be either due to an episodic outflow or
that it traces the dust condensation rim in the outflow. In both cases, we
think that the outer ring is linked with the outflow and probably with the
hourglass structure that we see at very large scales with ALMA. To distinguish
between these two interpretation we need to re-image this target with the VLTI
to look for changes from epoch to epoch.

##  5 An imaging Large Programme: INSPIRING

These first two images, of IRAS08544-4431 and HD101584, show very different
and complex morphologies that can be seen in the infrared at milliarcsecond
scales as it was also indicated by the VLTI/PIONIER snapshot survey[19]. To
build a solid understanding of the very close environment of these binaries,
we need to know wether those images are representative of this class of
objects. An imaging survey of these systems is perfect to provide such an
understanding. We have therefore initiated a Large Programme (PI: Kluska) that
we named INterferometric Survey of Post-agb bInaries with their RING
(INSPIRING) on the Very Large Telescope Interferometer (VLTI) using the
PIONIER and GRAVITY instruments. The main goal of this Large Programme is to
image eleven targets to see if the picture we have from the images of
IRAS08544-4431 and HD101584 is applicable to other such systems. The
scientific goals are the following:

  * •

reveal the 3-dimensional morphology of the disk inner rim and constrain
disk/binary interactions

  * •

obtain the absolute distances, luminosities, and masses of the systems by
combining the radial velocity orbit with the spatially resolved angular
position of the companion

  * •

look for direct evidence for accretion from the circumbinary disk

  * •

confirm the presence of a circum-secondary accretion disk in the sample

  * •

study the origin of extended flux component

  * •

reveal a possible disk evolutionary sequence

Those goals will be achieved by performing image reconstruction of the
continuum, locating the line absorption/emission from GRAVITY relative to the
continuum derived from PIONIER and modelling the dataset together with the
photometry with state-of-the-art 3-dimensional radiative transfer models[38].

For the image reconstruction part, the procedure is the following:

  * •

Perform an image reconstruction by subtracting a single star and fitting the
chromatic parameters ($f^{*}_{0}$ and $d_{\mathrm{env}}$) in the process. In
the case of a complete different morphology, perform classical image
reconstruction.

  * •

Perform the bootstrapping method to estimate the areas in the image where the
emission is significant to be taken into account

  * •

(optional) in case of detection of a second point source, perform another
image reconstruction by subtracting the second point source as well. This step
can be helped by model fitting to determine the flux and position of the
secondary.

  * •

Analyse the image geometry (e.g. is the ring centred on the primary?
secondary? Can we deduce the centre of mass?)

  * •

Use the obtained image to guide geometrical model fitting to the data and
derive sizes, temperatures and flux ratios between the different components.

Such an approach should be sufficient to extract the most important features
from the interferometric dataset and provide a direction for the
interpretation and further modelling using more physical models. A final
consistency check should be done by generating synthetic dataset with noise
from the final model and performed image reconstruction using the same hyper-
parameters and compare to the original image.

Additionally, more advanced image reconstruction features can be tested on
such a rich dataset such as data weighting (e.g., uniform or Briggs) that is
used in radio interferometry. Moreover, current image reconstruction
algorithms use generic regularisation functions that are not representing our
prior knowledge of the targets. We, therefore, have developed a new image
reconstruction algorithm based on deep learning techniques to learn a
regularisation function from images generated by astrophysical models. This is
presented in this volume of conference proceeding by Claes et al.

##  6 Conclusions

Infrared interferometry brings important information to decipher binary
evolution and the interaction of the binary with its close environment. While
first observations and snapshot surveys suggested a high complexity of the
near-infrared emission geometry, only the first two images, on IRAS08544-4431
and HD101584, entirely revealed their morphology thanks to advances in image
reconstruction techniques. The SPARCO technique, originally developed for
protoplanetary disks around young stars, showed to be very efficient to reveal
the circumbinary environment of these targets. Moreover, the bootstrapping
technique to estimate the significant emission in the image was proven to be
useful to guide the astrophysical interpretation and modelling of the images.
We also show the complementarity between ALMA and the VLTI. The image of
IRAS08544-4431 has particularly shaped our knowledge of these systems, by
spatially resolving the different building blocks and allowing to analyse them
one with respect to the other. An imaging Large Programme is on-going to test
this archetype and look for diversity in these intriguing and important to
understand evolved binaries. VLTI and CHARA are perfect infrared
interferometric observatories to images these systems and progress on our
knowledge.

###### Acknowledgements.

J.K. acknowledges support from FWO under the senior postdoctoral fellowship
(1281121N). J.K. and H.V.W. acknowledge support from the research council of
the KU Leuven under grant number C14/17/082. J.A. and V.B. acknowledge support
from the Spanish MICINN under project grants AYA2016-789994-P (AxiN) and
PID2019-105203GB-C21 (EVENTs / WEBULAE WEB).

## References

  * [1] Van Winckel, H., “Post-AGB Stars,” Annual Review of Astronomy and Astrophysics 41, 391–427 (Jan. 2003).
  * [2] Gezer, I., Van Winckel, H., Bozkurt, Z., De Smedt, K., Kamath, D., Hillen, M., and Manick, R., “The WISE view of RV Tauri stars,” MNRAS 453, 133–146 (Oct. 2015).
  * [3] Oomen, G.-M., Van Winckel, H., Pols, O., Nelemans, G., Escorza, A., Manick, R., Kamath, D., and Waelkens, C., “Orbital properties of binary post-AGB stars,” A&A 620, A85 (Dec. 2018).
  * [4] Nie, J. D., Wood, P. R., and Nicholls, C. P., “Predicting the fate of binary red giants using the observed sequence E star population: binary planetary nebula nuclei and post-RGB stars,” MNRAS 423, 2764–2780 (July 2012).
  * [5] Van Winckel, H., “Binary post-AGB stars as tracers of stellar evolution,” arXiv e-prints , arXiv:1809.00871 (Sep 2018).
  * [6] Bujarrabal, V., Alcolea, J., Van Winckel, H., Santander- García, M., and Castro-Carrizo, A., “Extended rotating disks around post-AGB stars,” A&A 557, A104 (Sept. 2013).
  * [7] Bujarrabal, V., Castro-Carrizo, A., Alcolea, J., and Van Winckel, H., “Detection of Keplerian dynamics in a disk around the post-AGB star AC Herculis,” A&A 575, L7 (Mar. 2015).
  * [8] Bujarrabal, V., Castro-Carrizo, A., Alcolea, J., Santander-García, M., Van Winckel, H., and Sánchez Contreras, C., “Further ALMA observations and detailed modeling of the Red Rectangle,” A&A 593, A92 (Sept. 2016).
  * [9] Bujarrabal, V., Castro-Carrizo, A., Alcolea, J., Van Winckel, H., Sánchez Contreras, C., and Santander- García, M., “A second post-AGB nebula that contains gas in rotation and in expansion: ALMA maps of IW Carinae,” A&A 597, L5 (Jan. 2017).
  * [10] Bujarrabal, V., Castro-Carrizo, A., Winckel, H. V., Alcolea, J., Contreras, C. S., Santander- García, M., and Hillen, M., “High-resolution observations of IRAS 08544-4431. Detection of a disk orbiting a post-AGB star and of a slow disk wind,” A&A 614, A58 (June 2018).
  * [11] Deroo, P., van Winckel, H., Min, M., Waters, L. B. F. M., Verhoelst, T., Jaffe, W., Morel, S., Paresce, F., Richichi, A., Stee, P., and Wittkowski, M., “Resolving the compact dusty discs around binary post-AGB stars using N-band interferometry,” A&A 450, 181–192 (Apr. 2006).
  * [12] Deroo, P., Acke, B., Verhoelst, T., Dominik, C., Tatulli, E., and van Winckel, H., “AMBER and MIDI interferometric observations of the post-AGB binary IRAS 08544-4431: the circumbinary disc resolved,” A&A 474, L45–L48 (Nov. 2007).
  * [13] Hillen, M., Verhoelst, T., Van Winckel, H., Chesneau, O., Hummel, C. A., Monnier, J. D., Farrington, C., Tycner, C., Mourard, D., ten Brummelaar, T., Banerjee, D. P. K., and Zavala, R. T., “An interferometric study of the post-AGB binary 89 Herculis. I. Spatially resolving the continuum circumstellar environment at optical and near-IR wavelengths with the VLTI, NPOI, IOTA, PTI, and the CHARA Array,” A&A 559, A111 (Nov. 2013).
  * [14] Hillen, M., Menu, J., Van Winckel, H., Min, M., Gielen, C., Wevers, T., Mulders, G. D., Regibo, S., and Verhoelst, T., “An interferometric study of the post-AGB binary 89 Herculis. II. Radiative transfer models of the circumbinary disk,” A&A 568, A12 (Aug. 2014).
  * [15] Hillen, M., de Vries, B. L., Menu, J., Van Winckel, H., Min, M., and Mulders, G. D., “The evolved circumbinary disk of AC Herculis: a radiative transfer, interferometric, and mineralogical study,” A&A 578, A40 (June 2015).
  * [16] Acke, B., Degroote, P., Lombaert, R., de Vries, B. L., Smolders, K., Verhoelst, T., Lagadec, E., Gielen, C., Van Winckel, H., and Waelkens, C., “Amorphous carbon in the disk around the post-AGB binary HR 4049\. Discerning dust species with featureless opacity curves,” A&A 551, A76 (Mar. 2013).
  * [17] Hillen, M., Van Winckel, H., Menu, J., Manick, R., Debosscher, J., Min, M., de Wit, W. J., Verhoelst, T., Kamath, D., and Waters, L. B. F. M., “A mid-IR interferometric survey with MIDI/VLTI: resolving the second- generation protoplanetary disks around post-AGB binaries,” A&A 599, A41 (Mar. 2017).
  * [18] Menu, J., van Boekel, R., Henning, T., Leinert, C., Waelkens, C., and Waters, L. B. F. M., “The structure of disks around intermediate-mass young stars from mid-infrared interferometry. Evidence for a population of group II disks with gaps,” A&A 581, A107 (Sept. 2015).
  * [19] Kluska, J., Van Winckel, H., Hillen, M., Berger, J. P., Kamath, D., Le Bouquin, J. B., and Min, M., “VLTI/PIONIER survey of disks around post-AGB binaries. Dust sublimation physics rules,” A&A 631, A108 (Nov 2019).
  * [20] Le Bouquin, J.-B., Berger, J.-P., Lazareff, B., Zins, G., Haguenauer, P., Jocou, L., Kern, P., Millan-Gabet, R., Traub, W., Absil, O., Augereau, J.-C., Benisty, M., Blind, N., Bonfils, X., Bourget, P., Delboulbe, A., Feautrier, P., Germain, M., Gitton, P., Gillier, D., Kiekebusch, M., Kluska, J., Knudstrup, J., Labeye, P., Lizon, J.-L., Monin, J.-L., Magnard, Y., Malbet, F., Maurel, D., Ménard, F., Micallef, M., Michaud, L., Montagnier, G., Morel, S., Moulin, T., Perraut, K., Popovic, D., Rabou, P., Rochat, S., Rojas, C., Roussel, F., Roux, A., Stadler, E., Stefl, S., Tatulli, E., and Ventura, N., “PIONIER: a 4-telescope visitor instrument at VLTI,” A&A 535, A67 (Nov. 2011).
  * [21] Monnier, J. D., Millan-Gabet, R., Billmeier, R., Akeson, R. L., Wallace, D., Berger, J.-P., Calvet, N., D’Alessio, P., Danchi, W. C., Hartmann, L., Hillenbrand, L. A., Kuchner, M., Rajagopal, J., Traub, W. A., Tuthill, P. G., Boden, A., Booth, A., Colavita, M., Gathright, J., Hrynevych, M., Le Mignant, D., Ligon, R., Neyman, C., Swain, M., Thompson, R., Vasisht, G., Wizinowich, P., Beichman, C., Beletic, J., Creech-Eakman, M., Koresko, C., Sargent, A., Shao, M., and van Belle, G., “The Near-Infrared Size-Luminosity Relations for Herbig Ae/Be Disks,” ApJ 624, 832–840 (May 2005).
  * [22] Lazareff, B., Berger, J.-P., Kluska, J., Le Bouquin, J.-B., Benisty, M., Malbet, F., Koen, C., Pinte, C., Thi, W.-F., Absil, O., Baron, F., Delboulbé, A., Duvert, G., Isella, A., Jocou, L., Juhasz, A., Kraus, S., Lachaume, R., Ménard, F., Millan-Gabet, R., Monnier, J. D., Moulin, T., Perraut, K., Rochat, S., Soulez, F., Tallon, M., Thiébaut, E., Traub, W., and Zins, G., “Structure of Herbig AeBe disks at the milliarcsecond scale . A statistical survey in the H band using PIONIER-VLTI,” A&A 599, A85 (Mar. 2017).
  * [23] Hillen, M., Kluska, J., Le Bouquin, J.-B., Van Winckel, H., Berger, J.-P., Kamath, D., and Bujarrabal, V., “Imaging the dust sublimation front of a circumbinary disk,” A&A 588, L1 (Apr. 2016).
  * [24] Kluska, J., Hillen, M., Van Winckel, H., Manick, R., Min, M., Regibo, S., and Royer, P., “The perturbed sublimation rim of the dust disk around the post-AGB binary IRAS08544-4431,” A&A 616, A153 (Aug. 2018).
  * [25] Kluska, J., Malbet, F., Berger, J.-P., Baron, F., Lazareff, B., Le Bouquin, J.-B., Monnier, J. D., Soulez, F., and Thiébaut, E., “SPARCO : a semi-parametric approach for image reconstruction of chromatic objects. Application to young stellar objects,” A&A 564, A80 (Apr. 2014).
  * [26] Kluska, J., Benisty, M., Soulez, F., Berger, J.-P., Le Bouquin, J.-B., Malbet, F., Lazareff, B., and Thiébaut, E., “A disk asymmetry in motion around the B[e] star MWC158,” A&A 591, A82 (June 2016).
  * [27] Kluska, J., Berger, J. P., Malbet, F., Lazareff, B., Benisty, M., Le Bouquin, J. B., Absil, O., Baron, F., Delboulbé, A., Duvert, G., Isella, A., Jocou, L., Juhasz, A., Kraus, S., Lachaume, R., Ménard, F., Millan-Gabet, R., Monnier, J. D., Moulin, T., Perraut, K., Rochat, S., Pinte, C., Soulez, F., Tallon, M., Thi, W. F., Thiébaut, E., Traub, W., and Zins, G., “A family portrait of disk inner rims around Herbig Ae/Be stars. Hunting for warps, rings, self shadowing, and misalignments in the inner astronomical units,” A&A 636, A116 (Apr. 2020).
  * [28] Gorlova, N., Van Winckel, H., Ikonnikova, N. P., Burlak, M. A., Komissarova, G. V., Jorissen, A., Gielen, C., Debosscher, J., and Degroote, P., “IRAS 19135+3937: an SRd variable as interacting binary surrounded by a circumbinary disc,” MNRAS 451, 2462–2478 (Aug. 2015).
  * [29] Bollen, D., Van Winckel, H., and Kamath, D., “Jet creation in post-AGB binaries: the circum-companion accretion disk around BD+46°442,” A&A 607, A60 (Nov. 2017).
  * [30] Bollen, D., Kamath, D., Van Winckel, H., and De Marco, O., “A spatio-kinematic model for jets in post-AGB stars,” A&A 631, A53 (Nov 2019).
  * [31] Bollen, D., Kamath, D., De Marco, O., Van Winckel, H., and Wardle, M., “Determining mass-accretion and jet mass-loss rates in post-asymptotic giant branch binary systems,” A&A 641, A175 (Sept. 2020).
  * [32] Kluska, J., Olofsson, H., Van Winckel, H., Khouri, T., Wittkowski, M., de Wit, W. J., Humphreys, E. M. L., Lindqvist, M., Maercker, M., Ramstedt, S., Tafoya, D., and Vlemmings, W. H. T., “VLTI/PIONIER reveals the close environment of the evolved system HD 101584,” A&A 642, A152 (Oct. 2020).
  * [33] Olofsson, H., Vlemmings, W. H. T., Maercker, M., Humphreys, E. M. L., Lindqvist, M., Nyman, L., and Ramstedt, S., “ALMA view of the circumstellar environment of the post-common-envelope-evolution binary system HD 101584,” A&A 576, L15 (Apr 2015).
  * [34] Olofsson, H., Vlemmings, W. H. T., Bergman, P., Humphreys, E. M. L., Lindqvist, M., Maercker, M., Nyman, L., Ramstedt, S., and Tafoya, D., “First detection of methanol towards a post-AGB object, HD 101584,” A&A 603, L2 (Jul 2017).
  * [35] Olofsson, H., Khouri, T., Maercker, M., Bergman, P., Doan, L., Tafoya, D., Vlemmings, W. H. T., Humphreys, E. M. L., Lindqvist, M., Nyman, L., and Ramstedt, S., “HD 101584: circumstellar characteristics and evolutionary status,” A&A 623, A153 (Mar 2019).
  * [36] Bakker, E. J., Lamers, H. J. G. L. M., Waters, L. B. F. M., and Waelkens, C., “The 218-day period of the peculiar late B-type star HD 101584.,” A&A 310, 861–871 (Jun 1996).
  * [37] Díaz, F., Hearnshaw, J., Rosenzweig, P., Guzman, E., Sivarani, T., and Parthasarathy, M., “Radial-Velocity Analysis of the Post-AGB Star, HD101584,” in [Binary Stars as Critical Tools & Tests in Contemporary Astrophysics ], Hartkopf, W. I., Harmanec, P., and Guinan, E. F., eds., IAU Symposium 240, 127 (Aug 2007).
  * [38] Min, M., Dullemond, C. P., Dominik, C., de Koter, A., and Hovenier, J. W., “Radiative transfer in very optically thick circumstellar disks,” A&A 497, 155–166 (Apr. 2009).