Topological polarization networking in uniaxial ferroelectrics

Y. Tikhonov^1,2, J. R. Maguire^3, C. J. McCluskey^3, J. P. V. McConville^3, A. Kumar^3, D. Meier^4, A. Razumnaya^2, J. M. Gregg^3, A. Gruverman^5, V. M. Vinokur^6,7,1* & I. Luk'yanchuk^1

The basic property of a ferroelectric is that however complex, the configuration of the polarization distribution has to maintain its internal electric neutrality. The neutrality implies that the polarization vector field is divergenceless, hence similar to the velocity field in incompressible liquids. This similarity manifests in, for example, the recently discovered nanoscale vortex-like topological excitations in ferroelectrics <cit.> and whirls in water streams.
The polarization vector is free to rotate
in the common pseudocubic perovskite ferroelectrics, which leads to a plethora of the observed topological excitations <cit.>.

In uniaxial ferroelectrics, in contrast, the only established topological excitations are domain walls (DWs), because the spontaneous polarization has only two stable orientations aligned with the anisotropy polar axis.
Moreover, the electroneutrality condition imposes strong additional restrictions.
Namely, the inclination of the DWs with respect to
the polar axis,
leads to head-to-head (H-H) or tail-to-tail (T-T)
collisions of polarization vectors belonging to the
oppositely-oriented domains, adjacent to the DWs from different sides, hence to the emergence of the
bound charges <cit.>. To avoid the charging, domain walls need to align
with the polarization direction.
As a result, domains usually spread across the entire sample from the lower to upper surface acquiring a shape of cylinders or stripes.
At the same time, the long but finite cigar-like domains have been
observed in several
ferroelectrics <cit.>.
Deviations from the electroneutrality at the loci of domains terminations were attributed to the screening by free charge carriers, compensating the bound charges <cit.>.

Here we demonstrate that uniaxial ferroelectrics can, in fact, harbor chargeless polarization  configurations of arbitrary complexity.
We establish that following the same fundamental laws of topological hydrodynamics  as topological excitations
in cubic perovskite ferroelectrics, the domains arrange themselves into an
ample-structured topological network of counterflowing mutually avoiding polarization fluxes.

The electroneutrality of ferroelectrics follows from the requirement of minimizing the electrostatic energy associated with the depolarisation fields.
To achieve it,
the polarization distribution tends to assume a structure with vanishing bound charges so that $\mathrm{div}\,\textbf{P} = 0$. The divergenceless nature of the polarization field is a
basic condition defining the physics of the spatially nonuniform ferroelectricity and resulting in the topological excitations. Therefore, by their topological characteristics, the polarization lines are identical to the streamlines of an ideal incompressible liquid, which enables us to employ the methods of the topological hydrodynamics<cit.>.

The insight is gained by comparing encountering polarization fluxes with the meeting of countercurrents in incompressible liquid that slightly deviate from their original paths to avoid the clash.
The counterflowing streams in the liquid correspond thus to the oppositely oriented polarization domains which bend to avoid the H-H and T-T collisions.
Then the neutrality is provided by the topological entwining of the encountered polarization line streams, separated by the single path-connected DW extended through the system. The energy cost related to the domain bending is far below the expense of the electrostatic energy associated with the onset of uncompensated bound charges.

The domain networking structure.

a, Lateral and vertical PFM data from the polar and non polar face of lead germanate, respectively.

b, An inner-corner cut of the rectangular uniaxial PGO sample.
The panel displays the simulated structure of interpenetrating polarization domains with the respective up- and down- polarization directions shown by white arrows. The magnitude of the z-component of the polarization, $P_z$, is quantified by the colour legend.
The upper and bottom near-surface regions host Landau fractal multidomain structure arising to compensate surface bound charges. The green contours depict the cross sections of the domain walls.

c, The charge map of the horizontal cross-section area shown in panel b. The upper plane displays the spontaneous polarization charge $\rho_z=-\partial P_z$, the up– and down- arrows in the upper plane showing the local directions of the $P_z$ polarization. Middle plane displays the charge $\rho_{xy}=-\partial_x P_x - \partial_y P_y$ induced by the transverse polarization whose local directions are shown by the in-plane arrows. The bottom plane displays the total resulting charge $\rho=-\mathrm{div}\,\mathbf{P}$ and demonstrates that $\rho_{xy}$ effectively screens $\rho_z$ so that the total charge, $\rho$, is vanishingly small in the vicinity of the inclined DWs.

d, Visualization of the complex network of polarization lines formed by the counter-flowing polarization fluxes corresponding to “up" and “down" domains.

To reveal the topology of the encountering polarization fluxes we explore the domain structure in the ferroelectric, lead germanate,
Pb$_5$Ge$_3$O$_{11}$ (PGO), which is a uniaxial ferroelectric featuring a three-fold polar axis along the [001] crystallographic axis<cit.>. This material exhibits a remnant polarization of P = 4.8 $\mu$C/cm$^2$, which develops through a second order phase transition at $T_c=177$ $^\circ$C where the space group symmetry changes from P (a paraelectric phase) to P$_3$ (a ferroelectric phase). A single crystal of PGO is grown by the Czochralski method with the seed crystal being pulled from the melt along the [001] direction. The samples are oriented and cut with the diamond saw to produce the parallel-plane plates of several square millimeters in area and with the thickness of about 0.5 mm. The sample plates are polished with the diamond paste to create surfaces with an optical quality finish. Cutting along the (001) and (100) planes yields the samples with the polar (z-cut) and non-polar (x-cut) surfaces.

The as-grown domain structure in the PGO is mapped using the piezoresponse force microscopy (PFM) <cit.> on the x- and z-cut surfaces as shown in Fig. 1a. The measurements were taken
using a Veeco Dimension 3100 AFM system with a Nanoscope IIIa controller and external lock-in amplifier. An AC bias of amplitude 5V and frequency 20kHz was applied to a Pt-Ir coated silicon AFM tip (Nanosensors, PPP-EFM) while in contact with the sample surface.
The out-of-plane polarization distribution was mapped by monitoring vertical displacement of the cantilever when measurements were performed on the polar face of the crystal, and the in-plane polarization map was obtained by monitoring the lateral PFM signal while scanning the non-polar surface.
The PFM image acquired on the non-polar surface shows a pattern of domains elongated along the [001] polar axis. The three-dimensional illustration of the as-grown domain arrangement displays the difference in characteristic domain size along the different
crystallographic directions, reflecting the uniaxial nature of ferroelectricity of the material. The PFM image obtained on the polar surface displays an as-grown pattern of irregular domains with antiparallel out-of-plane polarization directions.
In addition to the neutral 180$^\circ$ domain walls, there are also sections of the domain boundaries that correspond to the H-H and T-T arrangements.
This observation corroborates the previously reported domain walls features <cit.> by showing that PGO is an unusual example of a uniaxial proper ferroelectric that readily develops uncharged T-T and H-H domain walls.

The domain wall's topological structure.

a, The collision and subsequent splitting of the polarization streams in point $\mathbb{V}$  result in the scissoring pattern.

b, The top view of the domain wall.

c, The side view of the domain wall.

d, The isometric projection of the domain wall. The point V  corresponds a scissoring of the streams in a.
Panels b–d demonstrate that domain wall is a single path-connected surface separating the channels carrying the continuous oncoming streams of polarization fluxes piercing the system.
The color map visualizes the effective surface charge density, $\sigma$, quantified by the the colour legend. Importantly, the  value of $\sigma$ is substantially less than the polarization reversal $\Delta P_z \simeq \pm 10 \mu C/cm^2$ at the head-to-head and tail-to-tail domain walls because of the lateral screening effect.

e,  The enlarged vicinity of the saddle singular point $\mathbb{V}$. The pink disc stands for the piece of the plain tangent to the saddle surface in the point $\mathbb{V}$. The blue lines depict the polarization flow in the vicinity of the singularity.

To gain further insight into the internal domain structure in the PGO, we perform a quantitative study of the polarization distribution, using phase-field modeling of the bumping polarization fluxes,
see Methods for details.
The simulations are based on the minimization of the Ginzburg-Landau functional, F=$\int f dV$  with respect to the magnitude of the z-directed spontaneous polarization, $P$, playing the role of the order parameter describing the uniaxial ferroelectric material and to the electric potential, $\varphi$, determining the electric field $\mathbf{E}=-\nabla \varphi$.
The free energy density is

\begin{gather}
f=\frac{1}{2}A(T-T_c)P^{2}+\frac{b}{4}P^{4}+\frac{c}{6}P^{6}+\frac{G}{2}\left( \nabla P\right) ^{2}+P\partial _{z}\varphi -\frac{1}{2}\varepsilon _{0}\varepsilon _{i}\left( \nabla \varphi \right) ^{2}\,.
\label{GL}
\end{gather}

The total polarization including also a non-critical part induced by the depolarization field is given by $\mathbf{P}=P\mathbf{z}+\varepsilon_0 \varepsilon_i \mathbf{E}$, where we assume that the background dielectric constant $\varepsilon_i$ is isotropic. The numerical values of the coefficients in Eq. (<ref>) are obtained from fitting the experimental data reported in the preceding studies of the PGO <cit.>,
see Supplementary Fig. 1, $A=1.09\times 10^{7}$ C$^{-2}$m$^{2}$NK$^{-1}$,
$b\approx 0$; $c=3.3\times 10^{14}$ C$^{-6}$m$^{10}$N,
$G=8.26\times 10^{-10}$ C$^{-2}$m$^{4}$N, the vacuum dielectric permittivity
$\varepsilon _{0}=8.85\times 10^{-12}$ C$^{2}$m$^{-2}$N$^{-1}$, and the background non-critical part of the dielectric constant
The encountering of oncoming antiparallel domains is achieved by the special choice of boundary conditions for polarization at the top and at the bottom of the sample, see Methods.

Figure 1b  displays the simulated structure of the interpenetrating domains with the respective up- and down- polarization directions shown by white arrows.
The fair correspondence between the experimentally observed, Fig. 1a, and simulated, Fig. 1b, domain textures is clearly seen.

To display the spatial polarization distribution, we make an inner-corner cut of the simulated PGO sample, see Fig. 1b.
The structure reveals the self-similar configuration with branching domains.
The vertical cross-cuts demonstrate the inclined domain walls, the cuts crossing the meanders producing the famed “cigar-shaped" forms.
The horizontal cross-cut exposes the irregular pattern of cross-sections of the domains with the DWs highlighted by green lines. Another important feature is that the upper and bottom near-surface regions host the Landau fractal multidomain structure that forms to compensate surface bound charges<cit.>.

To summarise here, both experimental data and simulations establish the existence of the bent DWs, where, supposedly, substantial bound charges, destabilizing ferroelectricity should arise.
At the same time, the electric charge, $\rho=-\mathrm{div}\,\mathbf{P}$, appears to be vanishingly small in the vicinity of the DWs ridges, where it is supposed to be the maximal
(see the charge map of the horizontal cross-section area, magnified at the lower plane of the Fig. 1c).
The bound charge at the bent DWs is indeed much smaller than the charge  expected from the $z$-gradient of the spontaneous polarization, $\rho_z=-\partial_z
P_z$, presented at the upper plane of Fig. 1c; the up– and down- arrows in the upper plane show the directions of the $P_z$ polarization.
The solution of the puzzle is that the $\rho_z$ charge is nearly compensated by the shown in the middle plane charge $\rho_{xy}=-\partial_x P_x - \partial_y P_y$ induced by the transverse, $P_x$ and $P_y$, polarization components, that effectively screen the charges $\rho_z$; the in-plane arrows in the middle plane indicate the directions of the transverse polarization.
Hence, the total bound charge arising near the bent DWs, $\rho=\rho_{z}+\rho_{xy}$, comes out much smaller than $\rho_z$.
The more detailed distribution of the polarization fields and charges is presented in Supplementary Fig. 2.
This figure also reveals the softening of the H-H and T-T DWs resulting in their broadening  with respect to the vertical uncharged DWs which occurs due to the interplay between the polarization and electric field induced by the bound charges as predicted in <cit.>.

The obtained smallness of the bound charge at DWs, supports the basic concept that uniaxial ferroelectrics maintain the divergenceless character of the polarization. Indeed, Fig. 1d displays the polarization distribution as flux bundles visually identical to the flowing towards each other and clash-avoiding streams of the incompressible liquid. Accordingly, the polarization forms a network of the counterflowing slightly bent domains alternatively carrying “up" and “down" polarization fluxes.
It is the transverse component of the polarization that ensures that the domains with opposite flux avoid the mutual clash.

The counterflowing domains are separated by domain walls as shown in detail in Fig 2.
Figure 2a presents an essential building block of the domain network,
the two branching streams assuming the V-shape clash-avoiding scissoring  configuration.
The global image of the domain wall is shown in Fig. 2b and constitutes a path-connected topological manifold. Figure 2c shows the side view of the DWs manifold. Figure 2d is a tilted view clearly exposing the scissoring point $\mathbb{V}$ marked red.
This point realizes a saddle node singularity of the polarization vector field $\mathbf{P}$ (Fig. 2e).
The node is characterized by the singularity index calculated as the
sign of the Jacobian of the field $\mathbf{P}$ at point $\mathbb{V}$,  $N=\mathrm{sgn}\,\mathrm{det} \left\Vert \partial _{i}P_{j}\right\Vert$ <cit.>.
Accordingly, $N$ is equal to $+1$ when the outcoming in the $\pm$z-direction polarization fluxes of T-T encountering domains captivate the polarization flow converging to the $\mathbb{V}$-point from the x-y plane; this situation is shown in Fig. 2e. In the opposite case of the H-H domains encountering, $N=-1$.

Domain wall saddle points in lead germanate

a, Vertical PFM phase (left) of lead germanate at different depths into the crystal, obtained by PFM tomography, along with the domain wall trace inferred from the PFM (right). The scale bar is 2µm.

b, 3D domain wall reconstruction, obtained by triangulating the points identified as the the DW surface in b.

c, A zoomed in view of the domain wall saddle point highlighted in c. Arrows show the intersecting maxima and minima, and the color scale represents z height.

The singularity of the field at point $\mathbb{V}$ is a topological property. To see that, let $S^2$ be a closed sphere centered at $\mathbb{V}$ that is pierced by the polarization field $\mathbf{P}$. Consider now the map of $S^2$ onto the space of the tips of the normalized vector $\mathbf{n}=\mathbf{P}/P$ which is also a 2-dimensional sphere ${S^\prime}^2$. Thus, each point $\in S^2$ has an image on  ${S^\prime}^2$ corresponding to the vector $\mathbf{n}$ at this point. This correspondence is characterized by the degree of the map $S^2\overset{f}{\rightarrow} {S^\prime}^2$ <cit.>,

\begin{equation}
\mathrm{deg}\,f=\frac{1}{4\pi }\int_{S^2}\left( \mathbf{n}\cdot\left[ \partial _{\theta}\mathbf{n}\times \partial _{\varphi} \mathbf{n} \right] \right) {d\theta d\varphi}
\label{ind}
\end{equation}

where the integration is performed over the surface $S^2$, with $\varphi$ and $\theta$ being the spherical coordinates.
This quantity is an integer number giving the number of times that the image of sphere $S^2$ is wrapped over the sphere ${S^\prime}^2$, hence it is the topological characteristic of the point $\mathbb{V}$ that does not depend on the choice of the wrapping surface. Importantly, the degree of the map $S^2\overset{f}{\rightarrow} {S^\prime}^2$ for the sphere surrounding the point $\mathbb{V}$ is identical to the singularity index $N$ for this point,  $\mathrm{deg}\,f=N$. This fundamental statement establishes the deep relation between the topological properties of the point $\mathbb{V}$ and the differential properties of the polarization vector field, $\mathbf P$, in its vicinity. Therefore, the stability of the domain scissoring structure hosting the singular saddle point $\mathbb{V}$ is topologically protected.

It is clearly seen that in the close vicinity of the saddle point $\mathbb{V}$, the shape of the DW is the saddle surface i.e., is a surface that curves up and curves down in perpendicular directions. The color of the manifold, see Fig. 2d, represents the effective surface charge density at the manifold which remains indeed low.
The DW assumes the configuration that minimizes its energy under constraints imposed by the boundary conditions stemming from the domains' collisions. The shape of the DW
is similar to that of the soap film, which is attached to the outer fixed frame and shapes itself to minimize its surface energy <cit.>. The DW is thus an exemplary realization of a concept widely known in mathematics as the Lagrange minimal surface problem <cit.>.
Note that our multiconnected domain wall presents a more general case than the usual minimal surface
of the soap film: it accounts also for the minimization of the non-local energy of bending of the DW-generating counter-flowing polarization field lines.

One of the key characteristics of this DW surface is the presence of an array of the saddle
points seen clearly in Fig. 2b. Such saddle points, however cannot be detected using conventional 2D PFM domain
mapping. Fortunately, recent advances in tomographic
imaging  genuinely render high-resolution 3D images of ferroelectric microstructures <cit.> to access the 3D domain wall distribution in the PGO. Here we use the tomographic atomic force
microscopy (TAFM) <cit.>. The TAFM involves
rastering the tip of a scanning probe microscope, using significant compressive force, as high as
microNewtons, to sequentially mechanically remove thin layers from the specimen surface.
Imaging can be done either at the same time as the machining (and under the same conditions), or
in separate scans taken between machining runs, in which more conventional tip pressures can
be applied for imaging. We take the former approach. Machining and imaging in the same scanned
area, allows for acquiring a series of planar domain images at different depths into the sample to be collected,
as illustrated in Fig. 3a. Here the vertical PFM phase maps, taken at differing depths below a polished
PGO polar surface, are stacked alongside equivalent images in which the DWs are highlighted using image processing techniques in which maxima in colour gradients are

These images enable rendering of the DW surface in 3D (Fig. 3b).
Several important observations are made: firstly, the DW appears as a single
continuous meandering surface throughout most of the investigated area; secondly, the DW
surface contains multiple saddle points. The pink box in Fig. 3b highlights the most
obvious one found in this dataset. This region is isolated and magnified for clarity in Fig. 3c; additional saddle points are explicitly highlighted in the Supplementary Fig. 3.

Our results conclusively establish the existence of a new topological state of matter, corresponding to a counterflowing domain state where the DW separating clash-avoiding polarization counterflows is a unique multiconnected surface. This clash-avoiding structure characterized by a multitude of the saddle points at the DW surface allows to avoid formation of the uncompensated charges at the locations of the H-H and T-T polarization flow impacts
and maintains the electric neutrality of the interior of the ferroelectric material.  Our discovery of almost uncharged H-H and T-T regions in PGO is entirely consistent with the observation that the domain walls in PGO do not manifest enhanced conductivity, unlike systems in which bound charges emerging at DWs are screened by the intrinsic current-conducting carriers.
The reported results validate the introduced topological concept describing ferroelectrics in terms of incompressible fluids.

§ METHODS

§.§ Tomographic atomic force microscopy (TAFM).

The TAFM experiments (Fig. 3) were performed on the polar face of the same lead germanate crystal, using an Asylum Research MFP-3D Infinity AFM system with an internal lock-in amplifier. The TAFM technique employed involves simultaneous PFM imaging of the domain structure and AFM machining. A conductive diamond-coated silicon probe was scanned over the surface at a speed of 75 $\mu$ms$^{-1}$. The probe has a high stiffness constant of 80 Nm$^{-1}$, which, when combined with a high AFM deflection setpoint, produces a significant force on the sample surface. This resulted in the sequential removal of layers of material at a rate of approximately 0.2 $\mu$m per scan, up to a total depth of  4 $\mu$m, over a 20 $\mu$m $\times$ 20 $\mu$m region. Domain information was collected simultaneously via PFM (AC bias amplitude 2V and frequency  3.0 MHz) allowing for the domain evolution through the depth of the crystal to be mapped. The locations of the domain walls at various depths were extracted by taking the gradient of the PFM phase maps collected during TAFM. Using these domain wall locations, a 3D reconstruction of the domain wall surface (Fig 3a) was created using the “alphaShape” triangulation function, implemented in MATLAB. All PFM and TAFM data were collected at room temperature under ambient conditions.

Phase-field method.
To perform numerical calculations we used phase-field method implemented within FEniCS software package <cit.>. Nonlinear system of partial differential equations to solve are coming from the variation of free-energy functional (<ref>) with respect to polarization $P$:

\begin{equation}
-\gamma \frac{\partial P}{\partial t} = \frac{\delta F}{ \delta P},
\label{FunctionalVariation}
\end{equation}

where $\gamma$ is the time-scale parameter, which value is irrelevant for current case and is taken to be equal unity.

The nonlinear system of equations (<ref>) is closed by the linear system corresponding to electrostatic Poisson equation:

\begin{equation}
\varepsilon_{0} \varepsilon_{i} \nabla^{2} \varphi = \partial_z P,
\label{Poisson}
\end{equation}

where $\varphi$ is the electrostatic potential.

The computational domain is a rectangular box with dimensions $200 \times 200 \times 500$ nm, that was partitioned into tetrahedrons with the help of 3D mesh generator gmsh <cit.>.
At the top and bottom surfaces of the computational domain the Dirichlet boundary conditions are imposed, $P_{top} = P^{\prime}\mathrm{sin}(2\pi x / 100 + \pi) \mathrm{sin}(2\pi y / 100)$, $P_{bottom} = P^{\prime}\mathrm{sin}(2\pi x / 100) \mathrm{sin}(2\pi y / 100)$, where $P^{\prime}=5$ $\mu$C cm$^{-2}$, the value close to the equilibrium value of polarization throughout the volume of sample. The phase shift between polarization distributions at the top and bottom rectangular surfaces was properly selected to ensure encountering of the oppositely polarized domains in the bulk of the sample. In the $x$ and $y$ directions the boundary conditions on $P$ and $\varphi$ are set to be periodic.

Approximation of time derivative on the left hand side of variation (<ref>) is accomplished by BDF2 stepper with variable time step <cit.>.
At the first time step of the simulation, the values of the polarization in the computational domain are taken randomly from the $[-10^{-4},10^{-4}]$ $\mu$C cm$^{-2}$ interval.
Solution of nonlinear system of equations is performed with the Newton based nonlinear solver with line search and generalized minimal residual method with restart <cit.>. Linear system resulting from discretization of Poisson  equation (<ref>) is solved using the generalized minimal residual method with restart.

§ ACKNOWLEDGEMENTS
This work was supported by H2020 ITN-MANIC action (Y.T, J.M.G , I.L.). J.M.G., C.J.McC. and J.P.V.McC., acknowledge funding from the Engineering and Physical Sciences Research Council (grant number EP/P02453X/1 and PhD studentship support). J.M.G and J.R.M. acknowledge studentship support from the Northern Ireland Department for the Economy (DfE).  Research at
the University of Nebraska (A.G.) was supported by the National Science
Foundation (NSF) grant DMR-1709237. The work of V.M.V. was supported by Fulbright Foundation and by Terra Quantum AG.

Determining Ginzburg-Landau functional parameters from the experimental data <cit.>.

a, Fitting the experimental data (blue circles) of the polarization temperature dependence  reveals the critical exponent $1/4$ below $T_c$ which means that the coefficient $b$ in the functional (1) of the main text is anomalously small and the system is close to the tricritical point where order of the transition is changed from the 2nd to the 1st one <cit.>. The inset showing linear dependence of $P^4$ on $(T_c-T)$ offers a strong support to this fit.

b, Fitting the experimental temperature dependence of the dielectric constants, $\varepsilon_f$ and $\varepsilon_p$ of the ferroelectric and paraelectric phases, measured in the $z$-direction, reveals the Curie-like behaviour. The Curie constant at the paraelectric side  $C=10400$ K of the transition has the magnitude four times smaller then its magnitude at the ferroelectric side. This feature again reflects that the system is close to the tricritical point with $b\approx 0$. The upper inset, demonstrating the temperature dependence of the inverse dielectric constant, strongly confirms this relation as the angles between the tangents to $\varepsilon_f^{-1}(T)$ and $\varepsilon_i^{-1}(T)$ and temperature axis at $T=T_c$ differ by factor of 4. The downward bending of the $\varepsilon_f^{-1}(T)$ dependence at low temperatures is due to the background dielectric constant that is found to be equal $\varepsilon_i=21$. The bottom inset uses the magnified $\varepsilon$ scale to highlight the nearly temperature-independent behavior of the transversal dielectric constant appearing to be equal to $\varepsilon_i$, which is barely seen in the main plot. Identification of the Curie constant $C$ allows for determining parameter $A=(\varepsilon_0 C)^{-1}=1.09\times 10^{7}$ C$^{-2}$m$^{2}$NK$^{-1}$ in functional (1) of the main text, which, in turn, allows for finding parameter $c=3.3\times 10^{14}$ C$^{-6}$m$^{10}$ N from the fit of the polarization dependence of panel a. The coefficient $G=8.26\times 10^{-10}$ C$^{-2}$m$^{4}$N in the gradient term of functional (1), which is assumed to be isotropic, is estimated by taking the coherence length  $\xi_0=[G/(T_cA)]^{1/2}$ of the order the elementary cell size 0.4 nm.

Polarization and charge distribution .

a,  A sketch of the inner-corner cut of the rectangular uniaxial PGO sample reproducing Fig. 1a of the main text visualizing the simulated domain configuration.
The panel displays the structure of the interpenetrating polarization domains with the respective up- and down- polarization directions shown by white arrows. The magnitude of the z-component of the polarization, $P_z$, is quantified by the colour legend. The upper and bottom near-surface regions display a Landau fractal multidomain branching structure which arises to compensate the surface bound charges. The green contours in the horizontal plane of the corner cut, depict the cross sections of the domain walls. The AB and CD lines represent the along- and across-domain cuts along which the distribution of the polarization component $P_z$ is taken.

b, The sketch of the distribution of the spontaneous polarization along the AB and CD cuts. Importantly, the profile of the of the charged DW arising in the head-to-head (H-H) and tail-to-tail (T-T) polarization collisions along the AB is more broader than the corresponding profile uncharged domain wall along CD. The softening of the H-H and T-T domain wall occurs due to the interplay between the polarization and electric field induced by the bound charges as was predicted in <cit.>.

c, The 3D map of the polarization and bound charge distribution. The lines represent the flow of the polarization field and the spacial color map shows the corresponding distribution of the charges produced by the spontaneous polarization, $\rho_z=-\partial_z P_z$. The red triangles mark the direction of the polarization field.

d, The distribution of the polarization lines and the bound charge in the magnified  xy-plane that screens the spontaneous bound charges $\rho_z$.

e, The illustration of the screening effect. The spontaneous bound charges $\rho_z=-\partial_z P_z$ are almost completely balanced by bound charges $\rho_{xy}=-\partial_x P_x-\partial_y P_y$ produced by the transversal polarization field with components $P_x$ and $P_y$, so that the resulting charge $\rho$ nearly vanishes.

Domain wall saddle points.

a, PGO domain wall reconstruction. The color bar represents $z$-height below the surface, and the boxed regions indicate further saddle points.

b-c, Domain wall trace extracted from the PFM tomography data. b occurs towards the top surface of the sample, and c is obtained after some milling, at a greater depth. The boxed region shows where domain walls merge to produce a saddle point.

d, The saddle point corresponding to the box in b-c, as represented by the 3D domain wall recreation.

e-g, h-j and k-m  follow the same pattern, for different saddle points.

Throughout the figure, purple boxes correspond to saddle points visible from the top down (where domains merge with increasing depth), and orange boxes correspond to saddle points visible from the bottom up (domains are initially merged and then separate with increasing depth).

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

Nature Publishing Group

titleObservation of polar vortices in oxide superlattices

Taylor & Francis

journaltitleIntegrated Ferroelectrics
titleStability of vortex phases in ferroelectric easy-plane nano-cylinders

Nature Publishing Group

journaltitleNature Commun.
titleDiscovery of stable skyrmionic state in ferroelectric nanocomposites

Nature Publishing Group

titleObservation of room-temperature polar skyrmions

Nature Publishing Group

journaltitleNature Materials
titlePolar meron lattice in strained oxide ferroelectrics

Nature Publishing Group

journaltitleNature Commun.
titleHopfions emerge in ferroelectrics

titleDomains in ferroic crystals and thin films

journaltitlePhys. Rev.
titlePyroelectricity, internal domains, and interface charges in triglycine sulfate

American Institute of Physics

journaltitleJournal of Applied Physics
titleElectrostatic considerations in BaTiO3 domain formation during polarization reversal

Springer Science & Business Media

titleTopological methods in hydrodynamics

American Institute of Physics

journaltitleApplied Physics Letters
title5PbO$\cdot$3GeO$_2$ crystal; a new ferroelectric

The Physical Society of Japan

journaltitleJournal of the Physical Society of Japan
titleFerroelectricity in Pb$_5$Ge$_3$O$_{11}$

American Institute of Physics

journaltitleJourn. of Appl. Phys.
titleFerroelectric and optical properties of Pb$_5$Ge$_3$O$_{11}$ and its isomorphous compound Pb$_5$Ge$_2$SiO$_{11}$

Nature Publishing Group

journaltitleNature communications
titlePiezoresponse force microscopy and nanoferroic phenomena

Wiley Online Library

journaltitleAdvanced Functional Materials
titleObservation of unconventional dynamics of domain walls in uniaxial ferroelectric lead germanate

Taylor & Francis

titleDomain structure of lead germanate

Springer Science & Business Media

titleFerroelectric phenomena in crystals: physical foundations

journaltitleMaterials Science and Engineering: B
titleModeling of ferroelectric domains in thin films and superlattices

journaltitlePhys. Rev. Lett.
titleUniversal properties of ferroelectric domains

Springer Science & Business Media

titleModern geometry—methods and applications: Part II: The geometry and topology of manifolds

journaltitlearXiv preprint arXiv:1807.05200
titleSoap films with gravity and almost-minimal surfaces

booktitleMinimal Surfaces
titleMinimal surfaces

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleThickness scaling of ferroelectricity in BiFeO3 by tomographic atomic force microscopy

ACS Publications

journaltitleACS Applied Nano Materials
titleReal-time 3D imaging of nanoscale ferroelectric domain wall dynamics in lithium niobate single crystals under electric stimuli: implications for domain-wall-based nanoelectronic devices

titleAutomated Solution of Differential Equations by the Finite Element Method

journaltitleInt. J. Numer. Methods Eng.
titleGmsh: A three-dimensional finite element mesh generator with built-in pre- and post-processing facilities

journaltitleJ. Comput. Appl. Math.
titleAdaptive stiff solvers at low accuracy and complexity

titlePETSc Web page

Argonne National Laboratory

numberANL-95/11 - Revision 3.15
titlePETSc Users Manual
eading=none] to remove section "References"


# Disentangling Causal Effects from Sets of Interventions in the Presence of
Unobserved Confounders

Olivier Jeunen
Amazon
Edinburgh, UK &Ciarán M. Gilligan-Lee
Spotify & UCL
London, UK &Rishabh Mehrotra∗
Sharechat
London, UK &Mounia Lalmas
Spotify
London, UK Work done while author was at Spotify.

###### Abstract

The ability to answer causal questions is crucial in many domains, as causal
inference allows one to understand the impact of interventions. In many
applications, only a single intervention is possible at a given time. However,
in some important areas, multiple interventions are concurrently applied.
Disentangling the effects of single interventions from jointly applied
interventions is a challenging task—especially as simultaneously applied
interventions can interact. This problem is made harder still by unobserved
confounders, which influence both treatments and outcome. We address this
challenge by aiming to learn the effect of a single-intervention from both
observational data and sets of interventions. We prove that this is not
generally possible, but provide identification proofs demonstrating that it
can be achieved under non-linear continuous structural causal models with
additive, multivariate Gaussian noise—even when unobserved confounders are
present. Importantly, we show how to incorporate observed covariates and learn
heterogeneous treatment effects. Based on the identifiability proofs, we
provide an algorithm that learns the causal model parameters by pooling data
from different regimes and jointly maximizing the combined likelihood. The
effectiveness of our method is empirically demonstrated on both synthetic and
real-world data.

##  1 Introduction

The ability to answer causal questions is crucial in science, medicine,
economics, and beyond, see [Gilligan-Lee, 2020] for a high-level overview.
This is because causal inference allows one to understand the impact of
interventions.111Causal inference also allows one to ask and answer
_counterfactual_ questions, see [Perov et al., 2020] and [Vlontzos et al.,
2022]. In many applications, only a single intervention is possible at a given
time, or interventions are applied one after another in a sequential manner.
However, in some important areas, multiple interventions are concurrently
applied. For instance, in medicine, patients that possess many commodities may
have to be simultaneously treated with multiple prescriptions; in
computational advertising, people may be targeted by multiple concurrent
campaigns; and in dietetics, the nutritional content of meals can be
considered a joint intervention from which we wish to learn the effects of
individual nutritional components.

Disentangling the effects of single interventions from jointly applied
interventions is a challenging task—especially as simultaneously applied
interventions can interact, leading to consequences not seen when considering
single interventions separately. This problem is made harder still by the
possible presence of unobserved confounders, which influence both treatments
and outcome. This paper addresses this challenge, by aiming to learn the
effect of a single-intervention from both observational data and sets of
interventions. We prove that this is not generally possible, but provide
identification proofs demonstrating it can be achieved in certain classes of
non-linear continuous causal models with additive multivariate Gaussian
noise—even in the presence of unobserved confounders. This reasonably weak
additive noise assumption is prevalent in the causal inference and discovery
literature [Rolland et al., 2022, Saengkyongam and Silva, 2020, Kilbertus et
al., 2020]. Importantly, we show how to incorporate observed covariates, which
can be high-dimensional, and hence learn heterogeneous treatment effects for
single-interventions. Our main contributions are:

  1. 1.

A proof that without restrictions on the causal model, single-intervention
effects cannot be identified from observations and joint-interventions.
(§3.1,3.2)

  2. 2.

Proofs that single-interventions _can_ be identified from observations and
joint-interventions when the causal model belongs to certain (but not all)
classes of non-linear continuous structural causal models with additive,
multivariate Gaussian noise. (§3.2,3.3)

  3. 3.

An algorithm that learns the parameters of the proposed causal model and
disentangles single interventions from joint interventions. (§4)

  4. 4.

An empirical validation of our method on both synthetic and real data.222The
code to reproduce our results is available at
github.com/olivierjeunen/disentangling-neurips-2022. (§5)

##  2 Related Work

Disentangling multiple concurrent interventions: [Parbhoo et al., 2021] study
the question of disentangling multiple, simultaneously applied interventions
from observational data. They propose a specially designed neural network for
the problem and show good empirical performance on some datasets. However,
they do not address the formal identification problem, nor do they address
possible presence of unobserved confounders. By contrast our work derives the
conditions under which identifiability holds. We moreover propose an algorithm
that can disentangle multiple interventions even in the presence of unobserved
confounders—as long as both observational and interventional data is
available. Related work by [Parbhoo et al., 2020] investigated the
intervention-disentanglement problem from a reinforcement learning
perspective, where each intervention combination constitutes a different
action that a reinforcement learning agent can take. Unlike this approach, our
work explicitly focuses on modelling the interactions between interventions to
learn their individual effects. Closer to our work is [Saengkyongam and Silva,
2020], who investigate identifiability of joint effects from observations and
single-intervention data. They prove this is not generally possible, but
provide an identification proof for non-linear causal models with additive
Gaussian noise. Our work addresses a complementary question; we want to learn
the effect of a single-intervention from observational data and sets of
interventions. Additionally, another difference between our work and that of
[Saengkyongam and Silva, 2020] is that they do not consider identification of
individual-level causal effects given observed covariates. In a precursor to
the work by [Saengkyongam and Silva, 2020], [Nandy et al., 2017] developed a
method to estimate the effect of joint interventions from observational data
when the causal structure is unknown. This approach assumed linear causal
models with Gaussian noise, and only proved identifiability in this case under
a sparsity constraint. However, like [Saengkyongam and Silva, 2020], our
result does not need the linearity assumption, and no sparsity constraints are
required in our identification proof. Finally, others including [Schwab et
al., 2020, Egami and Imai, 2018, Lopez and Gutman, 2017, Ghassami et al.,
2021] explored how to estimate causal effects of a single categorical-, or
continuous-valued treatment, where different intervention values can produce
different outcomes. Unlike our work, they do not consider multiple concurrent
interventions that can interact.

Combining observations and interventions: [Bareinboim and Pearl, 2016] have
investigated non-parametric identifiability of causal effects using both
observational and interventional data, in a paradigm they call “data fusion.”
More general results were studied by [Lee et al., 2020], who provided
necessary and sufficient graphical conditions for identifying causal effects
from arbitrary combinations of observations and interventions. Recent work in
[Correa et al., 2021] explored identification of counterfactual—as opposed to
interventional—distributions from combinations of observational and
interventional data. Finally, [Ilse et al., 2021] investigated the most
efficient way to combine observational and interventional data to estimate
causal effects. They demonstrated they could significantly reduce the number
of interventional samples required to achieve a certain fit when adding
sufficient observational training samples. However, they only prove their
method theoretically in the linear-Gaussian case. In the non-linear case, they
parameterise their model using normalising flows and demonstrate their method
empirically. They only consider estimating single-interventions, and do not
deal with multiple, interacting interventions.

Additive noise models: While certain causal quantities may not be generally
identifiable from observational and interventional data, by imposing
restrictions on the structural functions underlying causal models, one can
obtain _semi-parametric_ identifiability results. One of the most common weak
restrictions used in the causal inference community are additive noise models
(ANMs), first studied in the context of causal discovery by [Hoyer et al.,
2009]—and still widely used today [Rolland et al., 2022]. ANMs limit the form
of the structural equations to be additive with respect to latent noise
variables—but allow nonlinear interactions between causes. [Janzing et al.,
2009] used ANMs to devise a method for inferring a latent confounder between
two observed variables. This is otherwise not possible without additional
assumptions on the underlying causal model. See [Lee and Spekkens, 2017], [Lee
et al., 2019], and [Dhir and Lee, 2020] for an extension of this approach
beyond ANMs. ANMs have also been employed by [Kilbertus et al., 2020] to
investigate the sensitivity of counterfactual notions of fairness to the
presence of unobserved confounding. Our work proves that in certain classes of
ANMs, the effect of a single-intervention can be identified from observational
data and sets of interventions—even in the presence of unobserved confounders.
Moreover, we show how to incorporate observed covariates in these ANMs to
learn the heterogeneous effects of single-interventions.

##  3 Model Identifiability

_Identifiability_ is a fundamental concept in parametric statistics that
relates to the quantities that can, or can not, be learned from data
[Rothenberg, 1971]. An estimand is said to be identifiable from data if it is
theoretically possible to learn this estimand, given infinite samples. If two
causal models coincide on said data then they must coincide on the value of
the estimand in question. Hence if one finds two causal models which agree on
said data, but disagree on the estimand, then the estimand is not identifiable
unless further restrictions are imposed. In this section, we provide
identification proofs for single-variable interventional effects from
observational data and joint interventions, for several model classes. Our
theoretical analysis provides insights into the fundamental limitations of
causal inference—and the assumptions that are required for identification.

#### Problem Definition.

We adopt the Structural Causal Model (SCM) framework as introduced by [Pearl,
2009]. An SCM $\mathcal{M}$ is defined by
$\langle\\{\bm{C},\bm{X},Y\\},\bm{U},\bm{f},\mathsf{P}_{U}\rangle$, where
$\\{\bm{C},\bm{X},Y\\}$ are endogenous variables separated into covariates
$\bm{C}$, treatments $\bm{X}$, and the outcome $Y$, $\bm{U}$ are exogenous
variables (possibly confounders), $\bm{f}$ are structural equations, and
$\mathsf{P}_{U}$ defines a joint probability distribution over the exogenous
variables.

The SCM $\mathcal{M}$ also induces a causal graph—where vertices represent
endogenous variables, and edges represent structural equations. Vertices with
outgoing edges to an endogenous variable $X_{i}$ are denoted as the parent set
of this variable, or $\text{PA}(X_{i})$. Typically, the observed covariates
$\bm{C}$ causally influence the treatments as well as the outcome, and are a
part of this set. Every endogenous variable $X_{i}$ (including $Y$) is then a
function of its parents in the graph $\text{PA}(X_{i})$ and a latent noise
term $U_{i}$, denoting the influence of factors external to the model:

$\displaystyle X_{i}$ $\displaystyle\coloneqq f_{i}(\text{PA}(X_{i}),U_{i}).$
(1)

In _Markovian_ SCMs, these latent noise terms are all mutually independent.
However, in general, distinct noise terms can be correlated according to some
global distribution $\mathsf{P}_{U}$. In this case, such correlation is due to
the presence of _unobserved confounders_.

An intervention on variable $X_{i}$ is denoted by $\text{do}(X_{i}=x_{i})$,
and it corresponds to replacing its structural equation with a constant, or
removing all incoming edges in the causal graph. The core question we wish to
answer in this work, is under which conditions the treatment effect of a
single intervention can be disentangled from joint interventions and
observational data. That is, given samples from the data regimes that induce

$\displaystyle\mathbb{E}[Y|X_{i}=x_{i},X_{j}=x_{j},C=c],\text{ and
}\mathbb{E}[Y|\text{do}(X_{i}=x_{i},X_{j}=x_{j}),C=c],$

when can we learn conditional average causal effects:

$\displaystyle\mathbb{E}[Y|\text{do}(X_{i}=x_{i}),X_{j}=x_{j},C=c],\text{ or
}\mathbb{E}[Y|X_{i}=x_{i},\text{do}(X_{j}=x_{j}),C=c]?$

In what follows, we first show that this quantity is not identifiable without
restrictions on the causal model—a proof by counterexample. We then go on to
prove, again by counterexample, that simply assuming ANMs without restrictions
on the structure of the causal graph—the core assumption made by [Saengkyongam
and Silva, 2020]—is not enough for this complementary research question.
Finally, we prove identifiability of the treatment effect for ANMs without
_causal_ interactions between treatments. Note that this latter case does not
mean treatments are independent: treatments can be influenced by observed
covariates and unobserved confounders, and their interactions on the outcome
are defined by the unrestricted function $f_{Y}$.

Table 1: SCMs for 3.1 $\mathcal{M}$ | $\mathcal{M}^{\prime}$
---|---
$X_{1}=U_{1}$ | $X_{1}=U_{1}$
$X_{2}=X_{1}U_{2}$ | $X_{2}=U_{2}$
$Y\>\>=X_{1}X_{2}U_{Y}$ | $Y\>\>=X_{1}X_{2}U_{Y}$

Table 2: Interventional distributions on $X_{1}$ under SCMs $\mathcal{M}$ and
$\mathcal{M}^{\prime}$.

$\mathsf{P}_{\mathcal{M}}(Y,X_{2}|\text{do}(X_{1}))$ | $Y=0$ | $Y=1$
---|---|---
$\text{do}(X_{1}=0)$ | $X_{2}=0$ | $1$ | $0$
$X_{2}=1$ | $0$ | $0$
$\text{do}(X_{1}=1)$ | $X_{2}=0$ | $1-p$ | $0$
$X_{2}=1$ | $0$ | $p$

$\mathsf{P}_{\mathcal{M}^{\prime}}(Y,X_{2}|\text{do}(X_{1}))$ | $Y=0$ | $Y=1$
---|---|---
$\text{do}(X_{1}=0)$ | $X_{2}=0$ | $1-p$ | $0$
$X_{2}=1$ | $p$ | $0$
$\text{do}(X_{1}=1)$ | $X_{2}=0$ | $1-p$ | $0$
$X_{2}=1$ | $0$ | $p$

###  3.1 Unidentifiability for general SCMs

For simplicity, but without loss of generality, we consider two treatment
variables $\\{X_{1},X_{2}\\}$ and no covariates. We will show that two
different SCMs $\mathcal{M}$ and $\mathcal{M}^{\prime}$ can yield identical
observational distributions, as well as joint interventional distributions.
They will also agree on the single-variable interventional distribution for
treatment $X_{2}$, but disagree on the single-variable effect of treatment
$X_{1}$. As such, given data from sets of interventions, this example shows
the treatment effect of single-variable interventions is not generally
identifiable. Our SCMs are defined in Table 1, where
$U_{1}=U_{2}=U_{Y}\sim\text{Bernoulli}(p)$. The observational, joint, and
single-variable interventional distributions on $X_{2}$ are identical under
$\mathcal{M}$ and $\mathcal{M^{\prime}}$; we defer them to the supplemental
material. The interventional distribution on $X_{1}$ differs for $\mathcal{M}$
and $\mathcal{M}^{\prime}$, as shown in Table 2. This counterexample shows
that single-variable interventional effects are not identifiable for general
SCMs, given observational data and joint interventions.

###  3.2 Unidentifiability for general ANMs

An Additive Noise Model is an SCM where the influence of the latent noise
variables is restricted to be additive in the structural equations. That is,
Eq. 1 is restricted to the form:

$\displaystyle X_{i}$ $\displaystyle\coloneqq f_{i}(\text{PA}(X_{i}))+U_{i}.$
(2)

Following [Saengkyongam and Silva, 2020], we additionally assume the noise
distribution to be a zero-centered Gaussian with an arbitrary covariance
matrix: $\mathsf{P}_{U}\sim\mathcal{N}(0,\Sigma)$. Table 3 defines two SCMs
that satisfy these assumptions. $\mathcal{M}$ and $\mathcal{M}^{\prime}$ yield
identical observational, joint interventional, and single-interventional
distributions on $X_{2}$. However—they disagree on the causal effect of
intervening on $X_{1}$. An intuitive underlying reason for this, is that we
_can_ identify the expression
$\mathbb{E}[X_{2}|X_{1}=x_{1}]=f_{2}(x_{1})+\mathbb{E}[U_{2}|X_{1}=x_{1}]$,
but have no tools to disentangle the effects coming from the structural
equation, $f_{2}(x_{1})$, from those stemming from the additive noise,
$\mathbb{E}[U_{2}|X_{1}=x_{1}]$. As a result,
$\mathbb{E}[X_{2}|\text{do}(X_{1}=x_{1})]=f_{2}(x_{1})$ is not identifiable,
proving that general ANMs with unrestricted causal structures are insufficient
for disentangling treatment effects. In this example, and in general SCMs of
this structure, joint interventions _can_ be disentangled for the
_consequence_ treatment $X_{2}$, but not for the _causing_ treatment $X_{1}$:

######  Theorem 1 (Identifiability of disentangled conditional average
treatment effects in additive noise models with a causal dependency between
treatments).

Let
$\mathcal{M}=\langle\\{\bm{C},\bm{X},Y\\},\bm{U},\bm{f},\mathsf{P}_{U}\rangle$
be an SCM, where

$\displaystyle X_{i}$ $\displaystyle=f_{i}(\bm{C})+U_{i},\qquad\quad
X_{j}=f_{j}(\bm{C},X_{i})+U_{j},\qquad\quad Y=f_{Y}(\bm{C},\bm{X})+U_{Y},$

and $\mathsf{P}_{U}\sim\mathcal{N}(0,\Sigma)$. The estimand
$\mathbb{E}[Y|\text{do}(X_{j}),C]$ is identifiable from the conjunction of two
data regimes:  (1) the observational distribution, and  (2) the joint
interventional distribution on $(X_{i},X_{j})$.

We prove this by showing that the structural equations are identifiable from
the joint interventional regime—whereas the necessary (co-)variances are
identifiable from observational data. Formal proofs are deferred to the
supplementary material.

Table 3: SCMs for 3.2, where $(U_{1},U_{2},U_{Y})_{\mathcal{M}}\sim\mathcal{N}(0,\Sigma)$ and $(U_{1},U_{2},U_{Y})_{\mathcal{M}^{\prime}}\sim\mathcal{N}(0,\Sigma^{\prime})$. |  $\mathcal{M}$ | $\mathcal{M}^{\prime}$
---|---
$X_{1}=U_{1}$ | $X_{1}=U_{1}$
$X_{2}=X_{1}+U_{2}$ | $X_{2}=2X_{1}+U_{2}$
$Y\>\>=X_{1}X_{2}+U_{Y}$ | $Y\>\>=X_{1}X_{2}+U_{Y}$
$\text{ where }\Sigma=\begin{bmatrix}1&1&1\\\ 1&1&1\\\ 1&1&1\\\
\end{bmatrix},\text{ and }\Sigma^{\prime}=\begin{bmatrix}1&0&1\\\ 0&0&0\\\
1&0&1\\\ \end{bmatrix}.$

###  3.3 Identifiability for Symmetric ANMs

The crux of the problem of non-identifiability in Section 3.2 comes from the
fact that treatment $X_{1}$ has a direct causal effect on treatment $X_{2}$.
In many realistic applications, this might never occur. Treatments will
decidedly be correlated, but this can be encoded in the SCM via either the
observed covariates or the unobserved confounders. Then still, as we have no
restrictions on the functional form of the structural equations with respect
to the treatments (i.e. $f_{Y}$), treatments can still yield highly nonlinear
interaction effects on the outcome. As such, we now focus on the case where
the SCM is devoid of causal links between treatments, but includes observed
covariates as well as unobserved confounders. In this case, the structural
equations take the following form:

$\displaystyle X_{i}=f_{i}(\bm{C})+U_{i},\qquad\qquad\qquad
Y=f_{Y}(\bm{C},\bm{X})+U_{Y}.$ (3)

Identifiability in this setting requires additional assumptions about the
correlation between the covariates and confounders. That is, without any
assumptions on the nature of the relation between $C$ and $U$, our estimand is
_not_ identifiable. We include a proof (by counterexample) of this result in
the supplemental material. Intuitively, the problem is that allowing
$\mathbb{E}[U_{Y}|C]$ to be non-zero impacts the identifiability of $f_{Y}$
under joint interventions—as $f_{Y}$ also depends on $C$—whilst this is
crucial to allow for disentanglement. As we have shown additional restrictions
on the class of SCMs are necessary, we now assume independence between $C$ and
$U$. Future work can explore whether this can be relaxed to allow for linear
relationships, among others.

Fig. 1 visualises the structure of the causal graph we assume, when we have
$K$ treatment variables. We call such causal graphs “symmetric”. SCMs with
this structure occur in many application areas from economics to online
advertising and beyond. In SCMs of this form, the causal effect of all single-
variable interventions can be disentangled from observational and joint
interventional data:

$X_{1}$$X_{2}$$\dotsc$$X_{K}$$U$$Y$$C$ Figure 1: DAG for a symm. SCM with
covariates and unobserved confounders.

######  Theorem 2 (Identifiability of disentangled conditional average
treatment effects in additive noise models with symmetric structure).

Let
$\mathcal{M}=\langle\\{\bm{C},\bm{X},Y\\},\bm{U},\bm{f},\mathsf{P}_{U}\rangle$
be an SCM, where

$\displaystyle X_{i}$ $\displaystyle=f_{i}(\bm{C})+U_{i},\quad\forall
i=1,\dotsc,K,$ $\displaystyle Y$ $\displaystyle=f_{Y}(\bm{C},\bm{X})+U_{Y},$

$C\perp\\!\\!\\!\\!\perp U$, and $\mathsf{P}_{U}\sim\mathcal{N}(0,\Sigma)$.
The estimand $\mathbb{E}[Y|\text{do}(X_{i}),C]$ is identifiable from the
conjunction of two data regimes:  (1) the observational distribution, and  (2)
any interventional distribution on a set of treatments
$\bm{X}_{\text{int}}\subseteq\bm{X}$ that holds $X_{i}$:
$X_{i}\in\bm{X}_{\text{int}}$.

$Y$$X_{1}$$X_{2}$$U$$\text{do}(X_{1}=x_{1})\qquad$$\qquad\text{do}(X_{2}=x_{2})$$C$
(a) Only $\mathbb{E}[Y|C=c,\text{do}(X_{2}=x_{2})]$ is generally identifiable.

$X_{1}$$X_{3}$$U$$X_{2}$$Y$$C$ (b) Only
$\mathbb{E}[Y|C=c,\text{do}(X_{3}=x_{3})]$ is generally identifiable.

$Y$$X_{1}$$X_{2}$$U$$\text{do}(X_{1}=x_{1})\qquad$$\qquad\text{do}(X_{2}=x_{2})$$C$
(c) All $\mathbb{E}[Y|C=c,\text{do}(X_{i}=x_{i})]$ are generally identifiable.

Figure 2: Causal DAGs illustrating under which conditions the single-variable
causal effect on the outcome $\mathbb{E}[Y|C=c,\text{do}(X_{i}=x_{i})]$ is
identifiable from observational and joint interventional data.

In this section, we have studied for several classes of SCMs whether causal
effects can be disentangled. Our results provide insights into the fundamental
limits of learning and inference, and help crystallise which assumptions are
necessary and sufficient to make disentanglement of effects from sets of
interventions feasible. Figure 2 visualises and summarises our key results
denoting _which_ single-variable causal effects can be disentangled in ANMs
with zero-mean Gaussian noise. In what follows, we present a learning
methodology and validate our theoretical identifiability results with
empirical observations on synthetic and real data.

##  4 Methodology

In this section, we introduce our methodology for estimating SCMs and
providing estimates for treatment effects under _any_ set of interventions.
Estimating an SCM from a combination of observational and interventional data
boils down to  (1) estimating the structural equations, and  (2) estimating
the noise distribution.  We extend the Expectation-Maximisation-style
iterative algorithm proposed by [Saengkyongam and Silva, 2020] to handle
observed covariates, and to _disentangle_ causal effects instead of
_combining_ them. The resulting method is not limited to learning single-
variable causal effects, as the results from [Saengkyongam and Silva, 2020]
allow us to extend these newly learned single-variable effects to sets of
interventions. As a result, our method is able to generalise to sets of
interventions that were never observed in the training data, with the only
restriction that every variable that makes up the set was part of _some_
intervention set in the training data.

Say we intervene on a subset of treatments
$\bm{X}_{\text{int}}\subseteq\bm{X}$, and
$\bm{X}_{\text{obs}}\equiv\bm{X}\setminus\bm{X}_{\text{int}}$. In general, we
can decompose a causal query with interventions on $\bm{X}_{\text{int}}$ as
follows:

$\mathbb{E}[Y|\bm{C};\text{do}(\bm{X}_{\text{int}});\bm{X}_{\text{obs}}]=f_{Y}(\bm{C};\bm{X})+\mathbb{E}[U_{Y}|\bm{X}_{\text{obs}}].$
(4)

Samples consist of observed values for all endogenous variables. For
convenience, we denote a sample by $\bm{x}$. Suppose we have data from $d$
different data regimes, corresponding to $d$ different sets of interventions
(the empty set $\emptyset$ corresponds to the observational regime). The full
dataset consists of samples and their corresponding interventions
$\mathcal{D}=\\{(\bm{X}_{\text{int}};\bm{x})\\}$.

Algorithm 1 SCM Estimation for Symm. ANMs

1: Input: Dataset $\mathcal{D}$

2: Output: Parameter estimates $\widehat{\theta}$, $\widehat{\Sigma}$

3: Initialise $\widehat{\theta}$ and $\widehat{\Sigma}$

4: while not converged do

5: // Solve for $\theta$ with fixed $\widehat{\Sigma}$

6: Optimise log-likelihood in Eq. 6

7: // Solve for $\Sigma$ with fixed $\widehat{\theta}$

8: Estimate $\widehat{\Sigma}$ from
$\widehat{U}=\bm{x}-\bm{f}(\bm{x};\widehat{\theta})$ return
$\widehat{\theta},\widehat{\Sigma}$

#### Estimating the structural equations.

We parameterise the structural equations with $\theta$, denoted as
$f(\cdot;\theta)$. The Gaussian likelihood with covariance matrix $\Sigma$ is
denoted as $\mathsf{P}_{U}(\cdot;\Sigma)$.

The likelihood for a single endogenous variable $x_{i}$ is defined as
$L(x_{i};\theta,\Sigma)=\mathsf{P}_{U}(x_{i}-f_{i}(\text{PA}(x_{i});\theta);\Sigma)$.
The likelihood for a sample $\bm{x}$ is defined as the product of the
likelihoods for every endogenous variable that was _not intervened on_ in that
sample:

$L(\bm{x};\bm{X_{\text{int}}},\theta,\Sigma)=\prod_{x_{i}\in\bm{X}\setminus\bm{X}_{\text{int}}}L(x_{i};\theta,\Sigma).$
(5)

Naturally, the log-likelihood of the dataset is then:

$\ell(\mathcal{D};\theta,\Sigma)=\sum_{(\bm{X}_{\text{int}},\bm{x})\in\mathcal{D}}\log
L(\bm{x};\bm{X_{\text{int}}},\theta,\Sigma).$ (6)

In principle, any iterative optimisation procedure can be used to maximise Eq.
6 when we fix the covariance $\Sigma$. Typically, an appropriate method is
chosen with respect to the model parameterisation.

#### Estimating the noise distribution.

Following [Saengkyongam and Silva, 2020], we note that the maximum likelihood
estimate for the covariance matrix of a multivariate normal can be directly
computed from the sample. This allows for efficient closed-form computation of
this step. Additionally, note that the assumption of zero-mean Gaussian noise
simplifies the proof and yields an analytical solution for this step, but this
is not a general limitation. Indeed—what matters is that we can estimate the
conditional noise $\mathbb{E}[U_{Y}|\bm{X}_{\text{obs}}]$, which can in
principle be estimated via a different model. Algorithm 1 shows the full
procedure to estimate the parameters of a symmetric ANM. At inference time,
Eq. 4 allows us to estimate the outcome under any—even unseen—set of
interventions.

##  5 Empirical Validation

We now empirically validate the effectiveness of our method in estimating SCMs
from observational and joint-interventional data, and assess the accuracy of
the inferred outcomes under varying sets of interventions. We do this both on
synthetic data and semi-synthetic data based on real data from a medical
study. The research questions we wish to answer from empirical results, are
the following:

RQ1: Are we able to accurately estimate the outcome under previously unseen
sets of interventions?
RQ2: Are we able to accurately estimate the structural equations and the noise
distribution?
RQ3: How does our method handle varying levels of unobserved confounding?

###  5.1 Synthetic experiment

First we adopt a simulation setup, which gives us the freedom to vary the
_true_ underlying SCM and observe performance differences among competing
methods. The structural equations $f_{i},f_{Y}$ in Eq. 3 are polynomials with
second-order interactions to illustrate the effectiveness of the learning
method even when treatment effects are highly non-linear, and the optimisation
surface is highly non-convex. We use the well-known Adam optimiser in our
experiments [Kingma and Ba, 2014]. As a baseline learning method we adopt a
regression model that estimates the outcome from its direct treatments using
pooled data from different regimes. Here, the presence of unobserved
confounders severely hinders the method to provide accurate estimates under
varying sets of interventions, further motivating the need for causal
models.333The code to reproduce our results is available at
github.com/olivierjeunen/disentangling-neurips-2022.

We let $|\bm{C}|=|\bm{X}|=4$, yielding
$\theta_{i}\in\mathbb{R}^{11}\quad\forall f_{i}$ , and
$\theta_{y}\in\mathbb{R}^{37}$ for the parameterisation of $f_{y}$. The
parameters for the structural equations are randomly sampled from a uniform
distribution over $[-2,+2]$. The covariance matrix
$\Sigma\in\mathbb{R}^{5\times 5}$ is uniformly sampled from a uniform
distribution over $[-1,+1]$, and then ensured to be positive semi-definite
through Ledoit-Wolf shrinkage [Ledoit and Wolf, 2004]. We repeat this process
10 times with varying random seeds and report confidence intervals over
obtained measurements.

We vary the size of the training sample $|\mathcal{D}|\in\\{2^{i}\forall
i=5,\ldots,13\\}$, where the set of intervened treatments
$\bm{X}_{\text{int}}$ is one of
$\\{\emptyset;\\{X_{0},X_{1}\\};\\{X_{1},X_{2}\\};\\{X_{2},X_{3}\\};\\}$. As
such, no single-variable interventions and the majority of the possible sets
of interventions are never observed. From this data, we estimate an SCM using
the procedure laid out in Algorithm 1.

Then, for every possible ($2^{|\bm{X}|}$) set of interventions, we simulate
$10\,000$ samples and estimate the outcome based on our estimated SCM using
Eq. 4. We report the Mean Absolute Error (MAE) between our estimated outcome
and the true outcome. Note that, as the true outcome is Gaussian, the optimal
estimate is the location of that Gaussian, which will have an expected
deviation of $\sqrt{\frac{2}{\pi}}\sigma$ where $\sigma$ is the standard
deviation on the noise parameter $U_{Y}$.

#### Estimating outcomes (RQ1).

Figure 3 visualises the results from the procedure laid out above, increasing
the size of the training sample over the x-axis and reporting the MAE on the
y-axis. Note that both axes are logarithmically scaled, and the shaded regions
indicate 95% confidence intervals. The plots clearly indicate that our method
is able to provide accurate and close to optimal estimates for the outcome,
under _any_ possible set of interventions. Indeed, the plots labelled
_observational_ , _do_ $(X_{0},X_{1})$, _do_ $(X_{1},X_{2})$ and _do_
$(X_{2},X_{3})$ show the accuracy of our method to estimate the outcome under
interventions that occurred in the training data—but the remaining twelve
plots relate to previously _unseen_ data regimes. While our results allow us
to _disentangle_ joint interventions, the results of [Saengkyongam and Silva,
2020] allow us to _combine_ interventions. As such, our method incorporates
(and subsumes) theirs, in order to generalise to arbitrary sets of
interventions. In contrast, the regression method that is oblivious to
confounders fails to accurately estimate the outcome under any data
regime—even those that are seen in the training data, or those where
confounders yield no influence on the outcome (i.e. all treatments are
intervened on). As such, the reported results validate that the proposed
method provides an unbiased and consistent estimator for the disentangled
causal effect, learned from sets of interventions in the presence of
unobserved confounders.

#### Estimating SCMs (RQ2).

We visualise the results with respect to the accuracy of the parameter
estimates in Figure 4. These measurements are obtained using the same
procedure and runs as for RQ1. Increasing the size of the training sample over
the x-axis, the leftmost plot shows the accuracy of the estimated parameters
for $f_{Y}$. The rightmost plot shows the accuracy of the estimated covariance
matrix for the multivariate normal that defines the noise distribution. As the
baseline regression method is oblivious to the noise distribution, it is not
included in the latter plot. Both axes are logarithmically scaled, and the
shaded regions indicate 95% confidence intervals. We observe that our method
is able to efficiently and effectively estimate the underlying causal
model—empirically validating the identifiability results presented in this
work.

Figure 3: Mean Absolute Error for estimating outcome under varying sets of
interventions, for varying sample sizes, comparing baseline regression with
our method. We empirically show our estimator is unbiased and consistent,
exhibiting optimal predictions even at modest training sample sizes. Figure 4:
Mean Absolute Error for estimating the parameters of the underlying SCM, for
varying sample sizes, comparing baseline regression with our method. We
empirically validate our estimator is unbiased and consistent, giving accurate
estimates even at modest training sample sizes.

###  5.2 Semi-synthetic experiment based on real medical data

We now show our method can accurately estimate unseen sets of interventions
(RQ1) and demonstrate it is consistent under varying levels of unobserved
confounding (RQ3) using a semi-synthetic setup on real-world data from the
International Stroke Trial database [Carolei, 1997]. This was a large,
randomized trial of up to 14 days of antithrombotic therapy after stroke
onset. There are two possible treatments: the aspirin allocation dosage
$\alpha_{a}$, and the heparin allocation dosage $\alpha_{h}$. The goal is to
understand the effects of these treatments on a composite outcome, a
continuous value in $[0,1]$ predicting the likelihood of patients’ recovery.

Figure 5: Mean Absolute Error for estimating the impact of interventions on
treatment $\alpha_{a}$, learned only from observational and joint
interventional data. We vary the magnitude of the covariance between
treatments and outcome to assess the effect of varying confounding.

We partially follow the semi-synthetic setup laid out in Appendix 3 of [Zhang
and Bareinboim, 2021]. Specifically, we adopt their probability table for the
joint observational distribution of the covariates: the gender, age and
conscious state of a patient. This table was computed from the dataset to
reflect a real-world observational distribution. While [Zhang and Bareinboim,
2021] deal with single discrete treatments, we extend to the continuous and
multi-treatment case where treatments can be interpreted as varying dosages of
aspirin and heparin. This procedure is described in detail in the
Supplementary Material.

We add zero-mean Gaussian noise to the treatments and the outcome, where we
randomly generate positive semi-definite covariance matrices with bounded non-
diagonal entries—varying the limit on the size of the covariances in order to
assess the effect of varying confounding on our method. We repeat this process
5 times, and sample $512$ observational and $512$ joint-interventional samples
(both $\alpha_{a}$ and $\alpha_{h}$ to learn the SCM from). We evaluate our
learning methods on $5\,000$ evaluation samples to predict the outcome under a
single-variable intervention on the aspirin dose $\alpha_{a}$.

Figure 5 shows results with respect to the accuracy of the learnt model for
making predictions on the outcome under varying single-variable interventions
on the aspirin dosage $\alpha_{a}$. Naturally, as the regression method has no
way of accounting for the confounding effect, it fails dramatically. Our
proposed method is able to consistently approximate the optimal
performance—which is attained by predicting the expectation of the outcome
distribution and limited by its variance. For increasing levels of
confounding, the gap between the optimal performing oracle and our learnt
model increases slightly. As such, this semi-synthetic experiment validates
the effectiveness of our proposed learning method, which can predict outcomes
under unseen sets of interventions (RQ1) and is robust to varying levels of
confounding (RQ3), for a use-case based on real-world medical data.

##  6 Conclusions & Future Work

In this work, we motivated the need for methods that can disentangle the
effects of single interventions from jointly applied interventions. As
multiple interventions can interact in possibly complex ways, this is a
challenging task; even more so in the presence of unobserved confounders.
First, we proved that such disentanglement is not possible in the general
setting, even when we restrict the influence of the unobserved confounders to
be additive in nature. By restricting the structure of the causal graph to be
symmetric we showed identifiability _can_ be acquired. Additionally, we showed
how to incorporate observed covariates, and have empirically demonstrated our
method.

Future work will tackle the case where the noise distribution is not
restricted to a zero-mean multivariate Gaussian. As mentioned in Section 4,
this assumption provides a closed-form expression for the conditional
expectation of the outcome noise given the observed treatments—but universal
function approximators could be used to obtain similar guarantees in more
general model classes. Investigating whether the assumption of causal
independence between the unobserved confounders and covariates can be relaxed,
to allow for linear relationships say, provides another avenue to explore.

## References

  * [Bareinboim and Pearl, 2016] Bareinboim, E. and Pearl, J. (2016).  Causal inference and the data-fusion problem.  Proceedings of the National Academy of Sciences, 113(27):7345–7352.
  * [Carolei, 1997] Carolei, A. e. a. (1997).  The international stroke trial (ist): a randomised trial of aspirin, subcutaneous heparin, both, or neither among 19435 patients with acute ischaemic stroke.  The Lancet, 349(9065):1569–1581.
  * [Correa et al., 2021] Correa, J. D., Lee, S., and Bareinboim, E. (2021).  Nested counterfactual identification from arbitrary surrogate experiments.  arXiv preprint arXiv:2107.03190.
  * [Dhir and Lee, 2020] Dhir, A. and Lee, C. M. (2020).  Integrating overlapping datasets using bivariate causal discovery.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 3781–3790.
  * [Egami and Imai, 2018] Egami, N. and Imai, K. (2018).  Causal interaction in factorial experiments: Application to conjoint analysis.  Journal of the American Statistical Association.
  * [Ghassami et al., 2021] Ghassami, A., Sani, N., Xu, Y., and Shpitser, I. (2021).  Multiply robust causal mediation analysis with continuous treatments.  arXiv preprint arXiv:2105.09254.
  * [Gilligan-Lee, 2020] Gilligan-Lee, C. (2020).  Causing trouble.  New Scientist, 246(3279):32–35.
  * [Hoyer et al., 2009] Hoyer, P., Janzing, D., Mooij, J., Peters, J., and Schölkopf, B. (2009).  Nonlinear causal discovery with additive noise models.  In Twenty-Second Annual Conference on Neural Information Processing Systems (NIPS 2008), pages 689–696. Curran.
  * [Ilse et al., 2021] Ilse, M., Forré, P., Welling, M., and Mooij, J. M. (2021).  Efficient causal inference from combined observational and interventional data through causal reductions.  arXiv preprint arXiv:2103.04786.
  * [Janzing et al., 2009] Janzing, D., Peters, J., Mooij, J., and Schölkopf, B. (2009).  Identifying confounders using additive noise models.  In 25th Conference on Uncertainty in Artificial Intelligence (UAI 2009), pages 249–257. AUAI Press.
  * [Kilbertus et al., 2020] Kilbertus, N., Ball, P. J., Kusner, M. J., Weller, A., and Silva, R. (2020).  The sensitivity of counterfactual fairness to unmeasured confounding.  In Uncertainty in Artificial Intelligence, pages 616–626. PMLR.
  * [Kingma and Ba, 2014] Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * [Ledoit and Wolf, 2004] Ledoit, O. and Wolf, M. (2004).  A well-conditioned estimator for large-dimensional covariance matrices.  Journal of Multivariate Analysis, 88(2):365–411.
  * [Lee et al., 2019] Lee, C. M., Hart, C., Richens, J. G., and Johri, S. (2019).  Leveraging directed causal discovery to detect latent common causes.  arXiv preprint arXiv:1910.10174.
  * [Lee and Spekkens, 2017] Lee, C. M. and Spekkens, R. W. (2017).  Causal inference via algebraic geometry: feasibility tests for functional causal structures with two binary observed variables.  Journal of Causal Inference, 5(2).
  * [Lee et al., 2020] Lee, S., Correa, J. D., and Bareinboim, E. (2020).  General identifiability with arbitrary surrogate experiments.  In Uncertainty in Artificial Intelligence, pages 389–398. PMLR.
  * [Lopez and Gutman, 2017] Lopez, M. J. and Gutman, R. (2017).  Estimation of causal effects with multiple treatments: a review and new ideas.  Statistical Science, pages 432–454.
  * [Nandy et al., 2017] Nandy, P., Maathuis, M. H., and Richardson, T. S. (2017).  Estimating the effect of joint interventions from observational data in sparse high-dimensional settings.  The Annals of Statistics, 45(2):647–674.
  * [Parbhoo et al., 2021] Parbhoo, S., Bauer, S., and Schwab, P. (2021).  Ncore: Neural counterfactual representation learning for combinations of treatments.  arXiv preprint arXiv:2103.11175.
  * [Parbhoo et al., 2020] Parbhoo, S., Wieser, M., Roth, V., and Doshi-Velez, F. (2020).  Transfer learning from well-curated to less-resourced populations with hiv.  In Machine Learning for Healthcare Conference, pages 589–609. PMLR.
  * [Pearl, 2009] Pearl, J. (2009).  Causality.  Cambridge university press.
  * [Perov et al., 2020] Perov, Y., Graham, L., Gourgoulias, K., Richens, J., Lee, C., Baker, A., and Johri, S. (2020).  Multiverse: causal reasoning using importance sampling in probabilistic programming.  In Symposium on advances in approximate bayesian inference, pages 1–36. PMLR.
  * [Rolland et al., 2022] Rolland, P., Cevher, V., Kleindessner, M., Russel, C., Schölkopf, B., Janzing, D., and Locatello, F. (2022).  Score matching enables causal discovery of nonlinear additive noise models.  arXiv preprint arXiv:2203.04413.
  * [Rothenberg, 1971] Rothenberg, T. J. (1971).  Identification in parametric models.  Econometrica, 39(3):577–591.
  * [Saengkyongam and Silva, 2020] Saengkyongam, S. and Silva, R. (2020).  Learning joint nonlinear effects from single-variable interventions in the presence of hidden confounders.  In Peters, J. and Sontag, D., editors, Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI), volume 124 of Proceedings of Machine Learning Research, pages 300–309. PMLR.
  * [Schwab et al., 2020] Schwab, P., Linhardt, L., Bauer, S., Buhmann, J. M., and Karlen, W. (2020).  Learning counterfactual representations for estimating individual dose-response curves.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 5612–5619.
  * [Vlontzos et al., 2022] Vlontzos, A., Kainz, B., and Gilligan-Lee, C. M. (2022).  Estimating categorical counterfactuals via deep twin networks.  In UAI 2022 Workshop on Causal Representation Learning.
  * [Zhang and Bareinboim, 2021] Zhang, J. and Bareinboim, E. (2021).  Bounding causal effects on continuous outcome.  In Proceedings of the 35nd AAAI Conference on Artificial Intelligence.

## Appendix

In what follows, we include complete and formal mathematical proofs for the
theorems presented in the paper, as well as more information about our
experimental setup for the semi-synthetic experiment performed to answer RQ3.

##  Appendix A Proof of Theorem 2: Identifiability in Symmetric Additive Noise
Models

This section provides a proof for identifiabiltiy of single-variable causal
effects in symmetric additive noise models (Theorem 2).

###### Proof.

Our causal estimand is the effect of intervening on $X_{i}$. For notational
convenience, we assume $k=i$, i.e. the intervention is on the last treatment.
From the observational data regime, we can trivially obtain the joint
$\mathsf{P}(\bm{C},X_{1},\dotsc,X_{k},Y)$. As such, we can condition on the
covariates and remaining treatment variables and then marginalise to obtain
$\mathbb{E}[Y|\text{do}(X_{k}=x_{k}),C=c]$. We can rewrite our causal estimand
as follows:

$\begin{gathered}\mathbb{E}[Y|C=c,X_{1}=x_{1},\dotsc,X_{k-1}=x_{k-1},\text{do}(X_{k}=x_{k})]=\\\
f_{Y}(c,x_{1},\dotsc,x_{k})+\mathbb{E}[U_{Y}|C=c,X_{1}=x_{1},\dotsc,X_{k-1}=x_{k-1}].\end{gathered}$
(7)

From the joint interventional data regime, we have access to the following
expectation:

$\begin{split}\mathbb{E}[Y|C=c,\text{do}(X_{1}=x_{1},\dotsc,X_{k}=x_{k})]=f_{Y}(c,x_{1},\dotsc,x_{k}).\end{split}$
(8)

Subtracting Eq. 8 from Eq. 7 shows that we only need to provide
identifiability for the conditional expectation on the outcome noise, given
the remaining treatment variables and the observed covariates:

$\begin{split}\mathbb{E}[U_{Y}|C=&c,X_{1}=x_{1},\dotsc,X_{k-1}=x_{k-1}]=\Sigma_{u_{y}}\Sigma^{-1}_{u_{x}}\bm{u}_{x},\\\
\;\text{where
}\Sigma_{u_{y}}&=\begin{bmatrix}\sigma_{Y1}&\dotsc&\sigma_{Y(k-1)}\end{bmatrix},\\\
\Sigma_{u_{x}}&=\begin{bmatrix}\sigma^{2}_{1}&\dotsc&\sigma_{1(k-1)}\\\
\vdots&\ddots&\vdots\\\
\sigma_{(k-1)1}&\dotsc&\sigma^{2}_{k-1}\end{bmatrix},\\\ \text{and
}\bm{u}_{x}&=\begin{bmatrix}x_{1}-f_{1}(c)&\dotsc&x_{k-1}-f_{k-1}(c)\end{bmatrix}^{\intercal}.\end{split}$
(9)

Here, $\sigma_{ij}$ denotes the covariance between noise variables $U_{i}$ and
$U_{j}$, and $\sigma_{Yi}$ denotes the covariance between the outcome $U_{Y}$
and $U_{i}$. There are two types of unidentified factors in these expressions:
the structural equations, that is, the $f_{i}$’s encapsulated in $\bm{u}_{x}$
from Eq. 9, and the parameters of the noise distribution, which are
encapsulated in $\Sigma_{u_{x}}$ and $\Sigma_{u_{y}}$. We tackle these in what
follows:

#### Identifying the structural equations.

As a direct result of our model definition, we obtain
$\mathbb{E}[X_{i}|C=c]=f_{i}(c)$ from the observational data regime. This
follows because in the DAG of Figure 1, $C$ and $U_{i}$ are independent for
all $i$. As such, the structural equations (and thus $\bm{u}_{x}$ from Eq. 9)
are identifiable.

#### Identifying the noise distribution.

For any pair of treatment variables, we can obtain
$\mathbb{E}[X_{i}|C=c,X_{j}=x_{j}]=f_{i}(c)+\mathbb{E}[U_{i}|X_{j}=x_{j}]$.
The latter term in this expression can then be rewritten as
$\mathbb{E}[U_{i}|U_{j}=x_{j}-f_{j}(c)]$, for fixed values of $x_{j}$ and $c$.
Because we have shown the structural equations ($f_{i},f_{j}$) to be
identifiable, we have that $\mathbb{E}[U_{i}|U_{j}]$ is identifiable, which
gives us the covariance $\sigma_{ij}$. Hence, the entirety of $\Sigma_{u_{x}}$
is identifiable. The same procedure can be used for every entry of
$\Sigma_{u_{y}}$, and the covariances $\sigma_{Yi}$ are identifiable as a
result.

It follows naturally that
$\mathbb{E}[Y|C=c,X_{1}=x_{1},\dotsc,X_{k-1}=x_{k-1},\text{do}(X_{k}=x_{k})]$
is identifiable. As we have data from the observational regime, any
marginalisation of this query is identifiable as well, concluding the proof. ∎

##  Appendix B Unidentifiability under Unconstrained SCMs

Tables 4 hold the full distributions for the counterexample in Section 3.1,
showing that single-variable interventional effects are unidentifiable from
observational and joint interventional data for unconstrained SCMs.

Table 4: Distributions under both SCMs $\mathcal{M}$ and
$\mathcal{M}^{\prime}$.

(a) Observational joint distribution. $\mathsf{P}(X_{1},X_{2},Y)$ | $Y=0$ | $Y=1$
---|---|---
$X_{1},X_{2}=0,0$ | $1-p$ | $0$
$X_{1},X_{2}=0,1$ | $0$ | $0$
$X_{1},X_{2}=1,0$ | $0$ | $0$
$X_{1},X_{2}=1,1$ | $0$ | $p$

(b) Joint interventional distribution. $\mathsf{P}(Y|\text{do}(X_{1},X_{2}))$ | $Y=0$ | $Y=1$
---|---|---
$\text{do}(X_{1}=0,X_{2}=0)$ | $1$ | $0$
$\text{do}(X_{1}=0,X_{2}=1)$ | $1$ | $0$
$\text{do}(X_{1}=1,X_{2}=0)$ | $1$ | $0$
$\text{do}(X_{1}=1,X_{2}=1)$ | $1-p$ | $p$

(c) Interventional distribution on $X_{2}$. $\mathsf{P}(Y,X_{1}|\text{do}(X_{2}))$ | $Y=0$ | $Y=1$
---|---|---
$\text{do}(X_{2}=0)$ | $X_{1}=0$ | $1-p$ | $0$
$X_{1}=1$ | $p$ | $0$
$\text{do}(X_{2}=1)$ | $X_{1}=0$ | $1-p$ | $0$
$X_{1}=1$ | $0$ | $p$

Table 5: SCMs for C, where $(U_{C},U_{1},U_{2},U_{Y})_{\mathcal{M}}\sim\mathcal{N}(0,\Sigma)$ and $(U_{C},U_{1},U_{2},U_{Y})_{\mathcal{M}^{\prime}}\sim\mathcal{N}(0,\Sigma^{\prime})$. |  $\mathcal{M}$ | $\mathcal{M}^{\prime}$
---|---
$C\>\>=U_{C}$ | $C\>\>=U_{C}$
$X_{1}=C+U_{1}$ | $X_{1}=2C+U_{1}$
$X_{2}=C+U_{2}$ | $X_{2}=2C+U_{2}$
$Y\>\>=CX_{1}X_{2}+U_{Y}$ | $Y\>\>=CX_{1}X_{2}+U_{Y}$
$\text{ where }\Sigma=\begin{bmatrix}1&1&1&1\\\ 1&1&1&1\\\ 1&1&1&1\\\
1&1&1&1\\\ \end{bmatrix},\Sigma^{\prime}=\begin{bmatrix}1&0&0&1\\\ 0&0&0&0\\\
0&0&0&0\\\ 1&0&0&1\\\ \end{bmatrix}.$

##  Appendix C Unidentifiability when $C$ is not independent of $U$

Table 5 provides the counterexample mentioned in Section 3.3 (just after
Equation 3) showing that single-variable interventional effects are
unidentifiable from observational and joint interventional data when the
covariates $C$ are not independent of the unobserved confounders $U$ in the
ANMs we consider.

Indeed, in model $\mathcal{M}$, $C$ is correlated with the latent
$U_{1},U_{2}$, telling us this model is non-Markovian: that there is an
unobserved confounder correlating $C$ and $U_{1},U_{2}$. $C$ is thus not
independent of this unobserved confounder. This is not the case in model
$\mathcal{M}^{\prime}$ Moreover, we can see that $\mathcal{M}$ and
$\mathcal{M}^{\prime}$ are identical under interchange of $X_{1}$ and $X_{2}$,
and have the same joint and marginal distributions, as well as the same joint
interventional distribution.

As the two SCMs coincide on observational and joint interventional
distributions, this proves that the confounding distribution is not
identifiable from these data regimes alone, and single-variable interventions
as well as the full SCM are unidentifiable. This demonstrates the need for the
additional restrictions between $C$ and $U$ in Theorem 3.2.

##  Appendix D Proof of Theorem 1: Identifiability with Causally Dependent
Treatments

This section provides a proof for identifiability of single-variable causal
effects when there is a causal dependency among treatments (Theorem 1).
Observational and joint interventional data are not sufficient in this case to
identify causal effects on _all_ treatments – but we can identify the causal
effect of intervening on the _consequence_ treatment instead of the _causing_
treatment.

###### Proof.

Our causal estimand is the effect of intervening on $X_{j}$. We can rewrite
our causal estimand as follows:

$\begin{split}\mathbb{E}[Y|X_{i}=x_{i},\text{do}(X_{j}=x_{j}),C=c]=f_{Y}(c,x_{i},x_{j})+\mathbb{E}[U_{Y}|C=c,X_{i}=x_{i}].\end{split}$
(10)

From the joint interventional data regime, we have access to the following
expectation:

$\begin{split}\mathbb{E}[Y|C=c,\text{do}(X_{i}=x_{i},X_{j}=x_{j})]=f_{Y}(c,x_{i},x_{j}).\end{split}$
(11)

Subtracting Eq. 11 from Eq. 10 shows that we only need to provide
identifiability for the conditional expectation on the outcome noise, given
the observed value for treatment $X_{i}$ and the observed covariates $\bm{C}$:

$\begin{split}\mathbb{E}[U_{Y}|C=c,X_{i}=x_{i}]=\mathbb{E}[U_{Y}|U_{i}=x_{i}-f_{i}(c)]=\frac{\sigma_{Yi}}{\sigma_{ii}}(x_{i}-f_{i}(c)).\end{split}$
(12)

Here, the first step comes from our SCM definition, and the second step comes
from the fact that we assume the noise distribution to be a zero-centered
multivariate Gaussian. As such, we need to identify the function $f_{i}$, the
variance on the noise variable $U_{i}$, and the covariance between $U_{i}$ and
$U_{Y}$. We obtain $\mathbb{E}[X_{i}|C=c]=f_{i}(c)$ directly from the
observational data regime. This makes the noise variable
$U_{i}=X_{i}-f_{i}(\bm{C})$ identifiable. Now, as a result, we can identify
its variance $\sigma_{ii}$ and covariance $\sigma_{Yi}$ from the observational
data regime, which concludes the proof. ∎

##  Appendix E Consistency under Varying Levels of Confounding (RQ3)

In order to validate the consistency of our learning method under varying
levels of unobserved confounding (RQ3), we adopt data from the International
Stroke Trial database [Carolei, 1997]. We partially follow the semi-synthetic
setup laid out in Appendix 3 of [Zhang and Bareinboim, 2021]. Specifically, we
adopt their probability table for the joint observational distribution of the
covariates: the gender, age and conscious state of a patient. This table was
computed from the dataset to reflect a real-world observational distribution.
They deal with discrete treatments—which we extend to the continuous case
where treatments can be interpreted as varying dosages of aspirin and heparin.
We model the structural equation on the outcome as:

$\begin{split}Y=~{}&0.1S-0.1A+0.25(C-1)+\alpha_{a}+0.75\alpha_{h}-3SA-0.1S\alpha_{a}\\\
-&0.3A\alpha_{a}+0.1S\alpha_{h}+0.2A\alpha_{h}+0.3C\alpha_{h}-0.45\alpha_{a}\alpha_{h}.\end{split}$
(13)

This loosely reflects the same intuitions as laid out in [Zhang and
Bareinboim, 2021], where $S$ is the gender, $A$ the age, $C$ the conscious
state, $\alpha_{a}$ the aspirin dose and $\alpha_{h}$ the heparin dose. Note
that this does not reflect correct medical knowledge or insights—the goal is
merely to have a polynomial with second-order interactions that we can learn
to model. As described in the main text, we add zero-mean Gaussian noise to
the treatments and the outcome, where we randomly generate positive semi-
definite covariance matrices with bounded non-diagonal entries—varying the
limit on the size of the covariances in order to assess the effect of varying
confounding on our method. We repeat this process 5 times, and sample $512$
observational and $512$ joint-interventional samples (both $\alpha_{a}$ and
$\alpha_{h}$) to learn the SCM from. We evaluate our learning methods on
$5\,000$ evaluation samples to predict the outcome under a single-variable
intervention on the aspirin dose $\alpha_{a}$.

All experiments ran in a notebook on a laptop.