# Realization of Causal Representation Learning to Adjust Confounding Bias in
Latent Space

Jia Li1, Xiang Li2, Xiaowei Jia3, Michael Steinbach1, Vipin Kumar1
University of Minnesota, 1 Computer Science and Engineering; 2 Bioproducts and
Biosystems Engineering.
3 University of Pittsburgh, Computer Science.
1 {jiaxx213, stei0062, kumar001}@umn.edu,2<EMAIL_ADDRESS>3
<EMAIL_ADDRESS>
###### Abstract.

Applying Deep Learning (DL) models to graphical causal learning has brought
outstanding effectiveness and efficiency but is still far from widespread use
in domain sciences. In research of EHR (Electronic Healthcare Records), we
realize that some confounding bias inherently exists in the causally formed
data, which DL cannot automatically adjust. Trace to the source is because the
Acyclic Causal Graph can be _Multi-Dimensional_ , so the bias and causal
learning happen in two subspaces, which makes it unobservable from the
learning process.

This paper initially raises the concept of Dimensionality for causal graphs.
In our case, the 3-Dimensional DAG (Directed Acyclic Graph) space is defined
by the axes of causal variables, the Absolute timeline, and Relative
timelines; This is also the essential difference between _Causality_ and
_Correlation_ problems.

We propose a novel new framework Causal Representation Learning (CRL), to
realize Graphical Causal Learning in latent space, which aims to provide
general solutions for 1) the inherent bias adjustment and 2) the DL causal
models generalization problem. We will also demonstrate the realization of CRL
with originally designed architecture and experimentally confirm its
feasibility.

Causal Representation, Deep Learning, Multi-Dimensional Causality, Higher-
Dimensional Feature, Inherent Confounding Bias

##  1\. Introduction

Graphical causal learning is an important topic for multiple research fields,
like meteorology, biology, medicine, epidemiology, social sciences, and
policy-making (wood2015lesson, ; vukovic2022causal, ; ombadi2020evaluation, ).
It aims to uncover the underlying causal relationships from observational
data, which is naturally generated during the causal process. Better extracted
causal information from observations can help us to disentangle questions such
as “whether and when the flood will happen”, “what medicines are most
effective against complications”, or “what long-term factors determine the
population growth in future”.

The study of causality originated from classical statistics, which has
developed for decades and contributed greatly to establishing knowledge in
many areas. Technological improvement in collecting high-quality data has
brought essential challenges to statistical causal analytics: learning
graphical causality from data to be Bayesian networks is NP-hard. Meanwhile,
machine learning techniques show great potential in causal learning;
specifically, it has become a burgeoning topic in the Artificial Intelligence
(AI) field as part of bedrock (li2020causal, ).

Conventionally, causal analytics begins with the model specification with
prior knowledge, assuming the underlying distributions, thus may lead to model
approximation errors (often called model selection bias in statistics). A
great advantage of such classical modeling is to keep it interpretable, thus
making it easy to leverage domain knowledge (scheines1997introduction, ). On
the other hand, Deep Learning (DL) applications bring outstanding
effectiveness to causal learning. They usually do not require model
specifications but reduce modeling approximation errors automatically by
optimizing neural network connections. Some recent works greatly enhanced
their model efficiency by approximating the discrete structure constraint as a
continuous optimization function (zheng2018dags, ; zheng2020causal, ), thus
bypassing the NP-hard combinational searching scheme. However, in the
meantime, the black-box nature of DL diminishes the modeling interpretability.

Figure 1. Existing Causal Learning Framework and three Standard Causal
Assumptions (referred to in Section 1 & 3.1).

We summarise the current causal learning framework as Fig 1, to describe the
processes between what we have known (i.e., Causal Knowledge) and what we can
see (i.e., Data Variables). The modeling methods discussed above provide
different options for people to realize Causal Discoveries or Reconstructions.
Here causal reconstruction refers to the applications of known relations onto
observational data, such as estimating causal effects between variables known
to be causally related, establishing a synthetic data generator to reflect
causal knowledge, and so on. From the modeling perspective, a causal relation
equals a correlation with a certain temporal order, which is indistinguishable
from the ones ordered in spatial or other ways; and, even further, the ones
without order - Since the learning process is throughout a sequence of records
whose order matters is not dissimilar to the one without. In other words, the
dependencies we learned from the data unnecessarily indicate causal relations
but are highly likely to be hybrid, even though the training data is in
chronological order. A typical hybrid scenario is with hidden confounders: for
a pair of variables, if their common causes are unobservable, the learned
pairwise causation between them can be biased - Such a bias can even “make up”
a pseudo correlation when the two are independent of each other.
Unfortunately, such a pseudo-relationship can only be identified manually,
given prior knowledge. In short, the interpretation process (based on domain
knowledge) is indispensable for generating causal relations, which are
commonly presented as Directed Acyclic Graphs (DAGs). Thus, it is easy to see
that the black-box nature of DL impedes its popularity in domain science
applications. The neural networks often provide no meaningful insights into
their learned dependencies, although with accurate outputs.

In Fig 1, we did not locate or mark Causal Inference, a classical concept from
statistics. While in our opinion, it can be considered as a combination of
some manual interpretations (e.g., determining directions of causation) and
causal reconstructions (e.g., causal effect estimation).

Recently, some efforts have been made toward interpretable DL models in causal
learning. The main principle is to encode the known graphical causal structure
among data variables by inferring (part of) the neural network architecture
(luo2020causal, ). However, this strategy can hardly achieve general success.
It is still very difficult to determine how well the knowledge has been
utilized during the networking optimization or why it has not been. Our
previous work (li2020teaching, ) has confirmed that networking optimization
cannot be automatically guided by causal knowledge in some practical
occasions, e.g., in EHR (Electronic Healthcare Records) data, some medical
effects can be “invisible” to DL models. It can be attributable to some
inherent biases existing between Knowledge and Observations, which brings
unadjustable confounding bias into DL.

This paper will initially discuss the existence of such a bias (in Section
3.1), named as _Causal Representation Bias ( CRB)_. Based on our analysis, the
fundamental reason is that the causal graph represented as DAG can be _Multi-
Dimensional_ due to multiple independent timelines (in Section 3.2). Briefly,
the confounding bias and our learning process can happen in distinct subspaces
(or planes) of the multi-dimensional DAG, which makes the bias “invisible” to
DL models. Furthermore, we will propose a new framework _Causal Representation
Learning ( CRL)_ to be a general solution for the inherent biases adjustment
(in Section 3.3); and also demonstrate a realization of CRL (in Section 4)
with a proposed novel architecture to extract representations of causal
Variables (in Section 4.1-4.2) and causal Effects(in Section 4.3).

Specifically, CRL aims to realize graphical model learning (including
reconstruction and discovery) in latent space. For one thing, it can transform
causal relations from different spaces into a common one (i.e., the latent
feature space) to perform adjustment; For another, it can explicitly preserve
the model interpretations in latent space, generally utilizable by DL. It is
worth mentioning that the term “CRL” was initially raised in another work
(scholkopf2021toward, ), where Scholkopf et al. focus on knowledge
generalization (making low-level knowledge reusable when learning high-level
knowledge) based on certain problems from computer vision field. While this
paper aims to broaden this concept as a general framework to provide solutions
for the causal DL model generalization problem.

The proposed architecture aims to realize the higher-dimensional
representation extraction (in contrast to the traditional lower-dim one),
which is imperative to provide a sufficient degree of freedom for CRL. In the
feasibility experiments (in Section 5), we perform the latent-space graphical
reconstruction and heuristic discovery on a practical causal learning problem,
successfully confirming the effectiveness of the proposed CRL realization
method.

In summary, the principal contributions of this paper:

  * •

initially raises Acyclic Causal Graph’s _Multi-Dimensionality_ , leading to
inherent confounding bias (named CRB) that is not automatically adjustable by
DL causal models. It also defines the essential difference between _Causality_
and _Correlation_ problems.

  * •

summarizes the current causal learning framework, analyzes its limitation, and
proposes a new one (CRL) to solve inherent biases and the DL model
generalization problem.

  * •

proposes the architecture of higher-dimensional representation extraction to
realize CRL and experimentally confirms its feasibility.

##  2\. Related Work

###  2.1. Representation Learning in Causality

The concept of representation learning was proposed about 100 years ago,
including the unsupervised principal component analysis (PCA) in 1901 and
supervised linear discriminant analysis (LDA) in 1936. In recent years, the
deep network architectures for representation learning have turned out to be
widely applied (zhong2016overview, ). Data representation describes how much
the explanatory factors of variation behind the data have been revealed
(bengio2012unsupervised, ), and thus plays an important role in determining
the success of machine learning methods. In the quest for applicable
Artificial Intelligence (AI), more powerful methods are motivated to be
designed in this direction.

Particularly, for machine learning of causality, the use of representation was
raised very recently in the last year (scholkopf2021toward, ), where Scholkopf
et al. proposed the term “causal representation learning” (CRL), but not mean
to concentrate on defining it as a general framework. In this paper, they
point out that “most work in causality starts from the premise that the causal
variables are given”, and describes three critical challenges we are facing,
which also constitute our demands for CRL: 1) The robustness of modeling is
still weak; 2) The learned knowledge lacks reusable mechanisms that allow it
to be generalized; 3) We do not have effective ways to interpret deep learning
results to form new knowledge.

From the probabilistic modeling perspective, the question of causal feature
learning (i.e., representation extraction) is to recover a set of latent
factors that describe the dependent relations (interpreted as causations) over
the observed data variables. For data variable $x$ with latent factor $h$,
their joint distribution can be decomposed as $P(x,h)=P(x|h)P(h)$, where
$P(h)$ is a _prior_ and $P(x|h)$ is a _likelihood_. Thus a general latent
factor model can be inherently interpreted as a latent cause model
(bengio2012unsupervised, ) directed from $h$ to $x$, denoted as $h\rightarrow
x$. This principle forms the basis of the proposed _Causal Representation
Learning_ to be reasonable.

###  2.2. Deep Learning in Causality

In recent few years, DL has been applied to realize highly efficient graphical
causal models. The key novelty is to turn the discrete possible-DAGs-space
searching problem into a smooth function optimization question (luo2020causal,
). Specifically, these methods quantify the Acyclic-ness of the graph’s
adjacency matrix to be a continuous constraint, which established an elegant
mathematical connection between the classical combinatorial algorithms and the
continuous optimization techniques of machine learning (luo2020causal, ). For
instance, Lachapelle et al. described how to convert the DAG learning into a
pure neural network learning one in 2019 (lachapelle2019gradient, ); and Zheng
et al. transformed it into a generalized nonparametric problem that requires
no modeling assumptions in 2020 (zheng2020causal, ).

Based on these realizations, people have made efforts toward the interpretable
DL causal models. One way is to figure out conditional independence by finding
zero derivatives of the graphical causal functions (lachapelle2019gradient, ).
Another way is to learn it from the weights of the neural networks
(zheng2020causal, ): among all Markov causal paths that lead the direction
from $X$ to $Y$, if at least one of them has weight $0$, then $Y$ does not
depend on $X$. Such a weights-detecting methodology is also called a Neural
Architecture Search (NAS).

As we have confirmed in our previous work (li2020teaching, ), DL may not able
to fully utilize the known causal structures in practice, with the inherent
obstacle of CRBs. Therefore, the network architecture or weights unnecessarily
provide desired interpretations. We propose to leave the interpretation task
to be still with conventional models where the knowledge originally come from,
and let DL concentrate on representation learning tasks - This strategy is
from our newly proposed CRL framework introduced in the next Section.

###  2.3. Existing Causal Models

Conventionally, the causal models are established on observed data variables.
The traditional statistical methods include constraint-based ones and score-
based ones. The constraint-based causal discovery, including PC and Fast
Causal Inference (FCI), typically detects causation by exploiting conditional
independences among data. PC is based on a no-hidden-confounder assumption
(equivalent to Causal Sufficiency), while FCI is not. However, there exists a
significant deficiency: they can necessarily but not sufficiently identify
causations, where the set of possible causal structures satisfying the same
conditional independence is called Markov Equivalence Class
(glymour2019review, ). The score-based methods usually belong to heuristics.
For example, Greedy Equivalence Search (GES) directs its search over the space
of equivalence classes using a properly defined score function under the no-
hidden-confounder assumption.

The graphical causal models are generated from machine learning, often
presented as Functional Causal Models (FCMs) or Structural Equation Models
(SEMs) (glymour2019review, ). They assume the value of $Y$ to be a
deterministic function of its direct causes $X$ with noise term $\epsilon$,
$Y=f(X,\epsilon;\theta)$ where $\theta$ indicates the parameters set involved
in $f$. They aim to decompose a causal DAG into a product of simpler factors,
where the separation is by conditional independence. It has been shown that
they can distinguish different DAGs from an equivalence class and are good at
leveraging prior knowledge.

##  3\. Causal Representation Bias

This section will initially raise three relevant terms: Acyclic Causal Graph’s
Dimensionality, Causal Representation Bias (CRB), and Causal Representation
Learning (CRL).

The multi-dimensionality of the acyclic causal graph is the underlying reason
for the inherent confounding bias (CRB). There exist multiple independent
timelines in the causal graph, including the _Absolute_ timeline where the
causal learning process happens, and the other _Relative_ timelines, from
where the possibly happened confounding bias is inherently unobservable by
causal learning. We will discuss an example in Section 3.2 to illustrate why a
causal graph must be in 3-dimensions to stay acyclic; and emphasize that
existence of _Relative_ timelines fundamentally differentiates the _Causality_
problems away from _Correlation_ ones.

There may exist multiple relative timelines independent of each other. Hence,
the simplest way to solve all possible inherent biases is to transform all
causal relations into a common feature space and adjust the confounding biases
(maybe) jointly, then the causal learning process can be properly performed in
this space. We call this novel framework Causal Representation Learning (CRL).

###  3.1. Inherent Existence

In statistical analytics, the specified models usually reflect the best
guesses from domain scientists about the underlying causality based on the
assumption that the variables in observational data can fully represent the
model’s variables, denoted as Causal Sufficiency. In other words, it is to
assume the consistency between Knowledge and Observations. Fig 1 illustrates
two types of its violations: One comes from Observational Limit (Blue-shaded),
often presenting as the confounding biases we can be aware of; Another is from
Incomplete Knowledge (Green-shaded), which may be reflected as unexplainable
biases or other phenomenons, potentially containing newly observed relations
to form new knowledge.

Since such a confounding bias comes from inappropriate model assumptions, can
we conclude that the advanced DL methods, without requiring model
specifications, can learn the underlying causality adequately by leveraging
prior knowledge? Or, despite interpretation, can they at least provide
sufficiently trustable predictions? Here, we emphasize that not all
confounding biases can be adjusted as modeling approximation errors. Some
inherent biases are “invisible” to the causal learning process and cannot be
eliminated by improving the models. However, the required knowledge for
adjusting these inherent biases may have been given, the operation of which is
out of the scope of the current framework. And that is the motivation for us
to reconsider it.

Additionally, unlike Causal Sufficiency, we are usually concerned little about
possible violations of the Causal Markov condition (e.g., wrongly inversed
causal direction) or Causal Faithfulness assumption (e.g., sampling bias of
selecting training data), empirically adjusted in causal learning or
interpretation processes.

###  3.2. Dimensionality of Acyclic Causal Graph

A regular causal relation is often represented as a Directed Acyclic Graph
(DAG), which emphasizes “acyclic” to guarantee that any transaction over the
variables has to be closed instead of an endless loop. To learn about such a
relationship in a regular data matrix with variable values as columns, the
learning process usually goes over the rows, i.e., different _Instances_.
Particularly, a temporal model (e.g., RNN model) requires properly ordered
rows representing the time series and, accordingly, to learn the variables’
value changing along the timeline. Here, an _Instance_ often refers to an
individual or object with a series of rows describing its temporal changing.
Fig 2 displays two causal relations, where the causal graphs over _Variables_
are given at the top, and accordingly, the temporal causal graph of an
individual _Instance_ is provided below.

Figure 2. Examples of causal graphs on _Variables_ , and correspondingly an
individual _Instance_ changes its value along the timeline to reflect such a
causal relation. (a) Left Side: No need for a timeline to stay Acyclic, equal
to a Correlation; (b) Right Side: the Cycle makes the timeline a necessary
dimension for constructing a DAG without cycles.

In Fig 2 (a), the temporal learning process going through the rows of all
_Instances_ in a sequence is equal to a learning process taking each row as an
individual _Instance_. In other words, the data rows do not necessarily
indicate a timeline but can be assigned any other meaning without changing the
model results. So this relation can be seen as equal to a regular
_Correlation_ without temporal learning. In Fig 2 (b), we can consider the
temporarily changing _Instance_ to be a patient at risk of diabetes. The
causal graph on _Variables_ contains cycles because LDL has two values, before
and after the patient takes Statin. The probability of getting a T2D diagnosis
has also been changed accordingly, forming another cycle. In short, it
requires one more dimension, the natural timeline, to construct an acyclic
graph. So the _Causality_ in (b) is distinguishable from (a).

Furthermore, Fig 2 (b) can also illustrate why there may exist _Causal
Representation Bias_ (CRB). In causal learning, the variable LDL may be able
to represent its two roles simultaneously, whose difference appears to be the
_medical effects_. However, once some confounding bias happens in the Statin-
LDL cycle, say, some extraneous reason makes the no-Statin-taken happen much
more times than taking-Statin (which is the usual case in reality), the
current LDL values cannot accurately represent its roles but provide biased
_medical effects_ in the causal learning process.

Figure 3. The 3-Dimensional DAG of Fig 2 (b), with one more variable BP. Blue
indicates observable in AT-Plane, while red means only observable in RT-Plane.
The LDL-Statin cycle is curved for better displaying overlapped edges.

Commonly, a temporal model goes through time series along the nature timeline
and learns the relationship among variables. This process can be seen as on a
plane, defined by the axis of all causal variables and the axis of an
_Absolute_ timeline, where the “absolute” indicates the plane to be fixed and
immovable for the current learning task and commonly with nature values. We
briefly name this plane as AT-Plane and assign the nature dates (in M/D/Y) to
be the timeline values in this case (denoted by $t$).

To construct a DAG for Fig 2 (b), we need to involve an additional axis of a
timeline, as shown in Fig 3, which describes how the _medical effects_
temporally develops. We name it _Relative_ timeline accordingly, and the plane
it defines together with the axis of variables is called RT-Plane. Notably,
such two timelines are independent; in our case, their axes can even be
orthogonal.

Specifically, in Fig 2 the two timelines of an instance are different. The one
in (a) is _Absolute_ one since all instances consistently change their values
along it, in other words, this is the timeline to describe _Variables_ ; While
the one in (b) is a _Relative_ timeline, which can only be used for an
_Instance_ that we are interested in, i.e., a patient possibly taking Statin
in few days. We intend to emphasize that a _Relative_ timeline is for data
instances rather than data variables since RT-Plane aims to describe the other
_Roles_ of the data variables, which cannot be observable in AT-Plane. For
example, in Fig 3, the role of LDL-after-taking-Stain can only be seen from
the perspective of RT-Plane but not AT-Plane. A pseudo-population has to be
concreted to make such a role “visible”, in other words, to create a dataset
to represent this role specifically. In this case, we must select all
instances throughout the data, conditioned as five-day-length segments of the
time series (from $\tau$ to $\tau+4$) during which the patient possibly takes
Statin.

Therefore, we conclude that learning on RT-Plane has to be _Asynchronous_ from
the causal learning on AT-Plane, because the dataset it requests is not ready
to use but instead initially “scattered” throughout the data. Now, suppose
confounding bias exists in RT-Plane, then a pre-processing is required to
_Adjust_ such a bias on the data to make sure the variables in AT-Plane can
properly represent the _Roles_ , and subsequently, let the target causal
learning be properly performed. Accordingly, we have the definition below.

_Definition 1. Causal Representation Bias (CRB) is Confounding Bias in
Relative timeline, whose adjustment has to be Asynchronous from the causal
learning process on the Absolute timeline. _

Here we only discussed the scenario of one _Relative_ timeline for a typical
medical effect. However, various medicines have different effects separately,
which may bring multiple RT-Planes on demand and give the variable multiple
independent roles. Particularly, no _Relative_ timeline makes _Causality_
equivalent to _Correlation_.

_Theorom 1. An Acyclic Causal Graph must have one Absolute timeline and can
have multiple independent Relative timelines that make it Multi-Dimensional. _

Here, we initially discuss the causal graph in 3D space but would be very
cautious about defining the term _Causal Graph’s Dimensionality_. Given the
background of EHR data analysis, it seems safe to assume RT-Plans to be
orthogonal with the AT-Plane, where the RT-Plans may or may not be mutually
parallel (e.g., intersections may indicate associated medical effects).
However, many scientific fields with causality concerns have numerous prior
knowledge, potentially bringing various possibilities and broader
perspectives. Such a multi-dimensionality may come from factors other than
timelines (e.g., some spatial factors), lead to more complex plane
intersections, or be in higher dimensions like 4D causal graphs or even
further. Moreover, the mathematical meaning of multi-dimensional acyclic
causal graphs also deserves a lot of exploits.

###  3.3. Proposed CRL Framework

Under the current causal learning framework (as shown in Fig 1), the existence
of CRB can easily be overlooked since it only considers the _Absolute_
timeline and focuses on the learning tasks processed along it. This problem is
particularly evident in Deep Learning (DL) models. Compared to traditional
statistical analytics, DL is more advanced in performance but tends to
downplay the interpretation of modeling due to its black-box nature. However,
being aware of CRB usually depends on how well the domain knowledge has been
applied. To be specific, the known causal relations are required to be
identifiable from the model and thereby can be applied to check possible
biases. Without any structural interpretation, even though a DL model has
occasionally “understood” the given causal knowledge, its generalization will
still be a problem.

Another trouble with causal deep learning is that the global optimization
strategy brings high efficiency but probably increases the overall biases
simultaneously. Unlike traditional heuristics, which can deal with the local
confounding bias in each step regardless of other variables, the global DL
models wrap up the optimizations of all (pairwise) causations in the graph in
one run, possibly letting the local confounding bias impact the global
performance.

In recent years, the demands for interpretable causal models significantly
boosted cross-disciplinary applications, including the numerous efforts to
make domain knowledge utilizable to DL models. However, these efforts are
often case-by-case for practical problems and can not be generalized easily.
For example, in the previous work (li2020teaching, ), we extract a linearly
computable estimator to partially adjust confounding bias in medical effects,
but with the premise that the known causal structure is luckily to be
derivable.

Figure 4. The Proposed Framework: Causal Representation Learning (CRL). Red
indicates representation feature space; green is observational data space; and
blue means interpretable knowledge space for conventional causal analytics.

Accordingly, we propose a new framework to realize the generalizable
adjustment for the inherent CRBs (as in Fig 4), named as _Causal
Representation Learning_ (CRL). The basic principle is to separate the
Adjustment and Learning into two individual processes, finished at two
different spaces (since for CRB adjustments, they have to be). This way can
make DL concentrate on the learning processes at latent space and let the
existing causal analytics or knowledge navigate model interpretation.
Specifically, we aim to realize an autoencoder to extract the Variable
representations individually; then among which the learning tasks (including
both causal discoveries and reconstructions) can be performed effectively. In
short, we establish Graphical Models in a common latent space to “downgrade” a
multi-dimensional causality problem to correlation in it.

On the other hand, such Latent Factor Graphical Models can give DL models
explicit structural interpretations in feature space, thus, can be widely used
to realize DL model generalization. From an applications perspective, this
advantage could be wider useful than CRB adjustment since it also benefits
correlations.

##  4\. Realization of Causal Representation

As mentioned in Section 2.1, a general latent factor model $P(x,h)$ can be
inherently interpreted as a latent cause model $P(x|h)$. Based on this
principle, we can expect reasonably decomposite pairwise causation from $x$ to
$y$ (denoted as $x\Rightarrow y$) as a causal chain $x\rightarrow h\rightarrow
v\rightarrow y$ where $v$ is the latent factor of $y$, and $x\rightarrow h$
relies on a _posterior_ $P(h|x)$. Such that the dependency between $x$ and $y$
is represented as the _posterior_ distribution $P(v|h)$. Suppose we have the
parametrized model $v=f(h)$, then the tuple $(v,h,f)$ can be called a
_Representation_ of $x\Rightarrow y$ in latent space, consisting of _Variable
representations_ $v$ and $h$, and _Effect representation_ $f$.

More importantly, by “stacking” or “splitting” these posterior distributions,
we can build up a _Directed Graphical Model_ in latent space, which is
flexible to be generalized or reused. For example, given causation
$x\Rightarrow y\Rightarrow z$, we are interested in the indirect causal effect
of $x$ on $z$, decomposed as $x\rightarrow h\rightarrow v\rightarrow
t\rightarrow z$, where $t$ is the latent factor of $z$, and we want to know
$P(t|h)=P(t|v)\big{|}_{v=f(h)}$. Since $f$ is known, we only need to learn a
conditional _Effect representation_ $g$, such that $t=g(v)\big{|}_{v=f(h)}$ by
stacking the unknown _posterior_ $P(t|v)$ onto the known _posterior_ $P(v|h)$.

In the following, we will start from the proposed autoencoder architecture for
higher-dimensional _Variable representation_ extraction (4.1) and the original
design of its critical layers (4.2); introduce how to model _Effect
representation_ in latent space and stack them to realize causal
reconstruction (4.3), and heuristic discovery algorithm accordingly (4.4).

###  4.1. Variable Representation Architecture

The regular autoencoder intends to deal with large-volume objects, like images
with numerous pixels, and aims to extract their lower-dimensional features,
which is more like a distillation process. However, our appeal is to
individually convert causal variables onto latent space, where their
dependencies (i.e., causal effects) can be properly modeled. Thus, the common
space should be capable of representing all variables, even though they are
mutually independent (the worst case of causal learning). Besides, it should
provide a sufficient degree of freedom for possibly built causal models, so we
must also consider the graphical complexity. For example, in our experimental
dataset, lengths of the $10$ variables are between $1\sim 5$ and added up to
be $32$ in total (refer to Table 2); their average node degree will be around
$2\sim 3$. Initially, we approximate that $64$-dimensional latent space is
about to satisfy the requirement, but it turns out to be highly redundant. It
is probably because variable attributes are strongly related, and we gradually
reduced the dimension to $16$. Besides, we envision that refinement or further
feature learning can be helpful to make it more efficient.

The biggest challenge to realizing such an autoencoder is extending the input
variable vector’s length (i.e., dimensionality), which should sufficiently
encode the associative information among its values. More importantly, this
process has to be _Invertible_ , to achieve high accuracy for the output
vector, unlike the regular autoencoders without concern of reconstructing each
pixel. We design a pair of layers to finish such two symmetric processes,
named _Encrypt_ and _Decrypt_ , respectively located at the Input and Output
of the autoencoder, as shown in Fig 5.

Figure 5. Proposed autoencoder architecture for an individual causal Variable
$x$ representation extraction.

The Encrypt works as a “feature amplifier” for extracting higher-order
features from the input vector, while Decrypt exactly reverses this process.
How much the feature order should reach is on practical demand, such as a
_double-wise_ amplifying to extract the 2-dimensional associations, or even
higher like _triple-wise_ to encode 3-dimensional features further. In this
work, we only demonstrate the double-wise extraction. To be specific, all
possible combinations of two values from this vector will be “encrypted” using
a set of constants named _Key_ , generated by Encoder. Then the _Key_ will be
copied by the Decoder to “decrypt” its output back to a pair of such values.
For a pair, one _Key_ means one time amplifying, so with multiple _Keys_ , we
can get a group of vectors differently amplified from the same two values;
then concatenate them to form a much longer new vector to be a double-wise
feature of the input pair of values. Because for a variable with length $n$,
its double-wise features can be arranged as a $(n-1)*(n-1)$ sized matrix, we
use a square to represent it in Fig 5 (unnecessarily to be a 2-dimensional
feature).

For example, Fig 5 contains four _Keys_ , which means extracting four times
amplified features. As a metaphor, in this architecture, the role of _Keys_ is
similar to the filters in a regular autoencoder.

Figure 6. Encrypt and Decrypt with _Key_ $\mathbf{\theta=(s,t)}$.

###  4.2. Variable Encrypt and Decrypt

Suppose the input variable denoted as
$\overrightarrow{X}=(x_{1},x_{2},\ldots,x_{n})$ has length $n$. The Encrypt is
defined as a function $f(x_{j};x_{i},\theta)$, a transformation of $x_{j}$
using $x_{i}$ as weights and use the _Key_ $\theta$ as parameters, where
$i\neq j$ and $i,j\in\\{1,\ldots,n\\}$.

Let $\theta=(s,t)$, s.t. $s(x)=w_{s}x$ and $t(x)=w_{t}x$, given $w_{s}$ and
$w_{t}$. We define that $f(x_{j})=x_{j}\otimes exp(s(x_{i}))+t(x_{i})$, where
$\otimes$ denotes the element-wise multiplication, and let $y_{j}=f(x_{j})$ to
be the corresponding output for $x_{j}$. Accordingly, the inversed function
$f^{-1}$ in Decrypt is going to be $(y_{j}-t(y_{i}))\otimes exp(-s(y_{i}))$.
Fig 6 illustrates the two symmetric processes. Because calculation of $f^{-1}$
does not involve inversed function $s^{-1}$ or $t^{-1}$, the two basic
transformation functions $s$ and $t$ can be defined very flexible to be non-
linear.

We assemble all $f$ functions as a set $\mathcal{F}(X;\Theta)$, where $\Theta$
is the set of all $\theta$. Then the two layers of Encrypt and Decrypt can be
denoted as $Y=\mathcal{F}(X;\Theta)$, and $X=\mathcal{F}^{-1}(Y;\Theta)$
respectively. Source codes for these two layers are free for downloading, and
the autoencoder realization is also included
111https://github.com/kflijia/bijective_crossing_functions.git.

We want to highlight that the work of Dinh, L., et al inspires the design
above (dinh2016density, ). One can perform higher-order associations by the
same principle, like the triple-wise feature extraction.

###  4.3. Effect Representation and Stacking

Fig 7 displays the architecture of learning the causal Effect representation
$h\rightarrow v$ for the causation $x\Rightarrow y$ using the $RNN$ model,
which is to construct the causal chain $x\rightarrow h\rightarrow v\rightarrow
y$. The feature vectors $h$ and $v$ have been initialized as the extracted
Variable representations individually, which have the same dimensionality
($h,v\in\mathcal{R}^{L}$), but hardly be considered as in a _common feature_
space. In this Effect learning process, $h$ and $v$ will keep updated during
optimizations to get “closer” to each other. It can be seen that both $h$ and
$v$ are shifting in $\mathcal{R}^{L}$ space toward finding the proper
locations minimizing their distance and at the same time letting the $RNN$
sufficiently model their dependence $P(v|h)$, which represents the estimated
causal Effect $h\rightarrow v$ in $\mathcal{R}^{L}$.

Figure 7. Architecture of learning the $RNN$ model in latent space to build
the causal Effect representation $h\rightarrow v$.

There are three optimization tasks with three objective functions
respectively, which should be sequentially updated in each iteration:

  1. (1)

Update _posterior_ $P(h|x)$ and _likelihood_ $P(x|h)$ to minimize the
reconstruction error of $x$.

  2. (2)

Update _posterior_ $P(h|x)$ and _Effect posterior_ $P(v|h)$ to minimize the
reconstruction error of $y$.

  3. (3)

Update _posterior_ $P(v|y)$ and _likelihood_ $P(y|v)$ to minimize the
reconstruction error of $y$

Figure 8. Example of stacking two Causal Representations to fully-connect
$x,y$ and $z$.

“Stacking” two Causal Representations means that given one of them has been
well established, we aim to build the second one from either end of the first
one to establish a fully-connected 3-nodes causal graph, which can eventually
shape as a colliding, a confounding, or a chain. This is essential to ask how
to stack two _Effect posteriors_ , e.g., stack $P(t|v)$ onto $P(v|h)$ to form
the chain $h\rightarrow v\rightarrow t$. There exist four possible occasions
by combining two binary options: is the _prior_ at the bottom a _cause_ or a
_result_? And how will the _prior_ at the top be?

Suppose $y$ is the data variable whose latent space is where the stacking will
happen, and we have $x$ and $z$ at each side of $y$ to be the connecting
candidates, presented as the representations $h$ at left and $t$ at right. As
shown in Fig 8, the _prior_ $P(v)$ at the bottom can be either a _cause_ of
any other latent factor (output toward the right) or the _result_ caused by
$h$ (from the input Effect $P(v|h)$). While for the top layer, Fig 8 displays
two possible roles of $t$: to be the new _cause_ that inputs Effect
$P(\alpha|t)$ (green colored), or the new _result_ that takes output Effect
$P(\beta|t)$ (red colored). Although with different directions, in
optimization of the top layer, $t$ always represents the fixed $prior$ to
support updating the _posterior_ , $P(\alpha|t)$ or $P(\beta|t)$.

For each occasion, multiple inputs and outputs exist to be the potential
entrance and exit, respectively, where the selection would be on demand during
causal reconstructions or discoveries. For example, let $\mapsto$ be the
notation that connects entrance and exit, then, $P(v|h)\mapsto y$ means output
$P(y|x)$, $P(\alpha|t)\mapsto y$ means output $P(y|z)$; $y\mapsto P(\beta|t)$
means output $P(z|y)$ with the ground truth $y$ input, while $P(v|h)\mapsto
P(\beta|t)$ means output $P(z|y)\big{|}_{y=f(x)}$ where $f$ is the estimated
causal effect $x\Rightarrow y$.

From a geometric perspective, the desired common feature space is further
specified every time we complete a stacking in $\mathcal{R}^{L}$. Suppose $L$
is large enough to provide a sufficient degree of freedom throughout the
causal graph. In that case, we can finally obtain a group of feature vectors
in $\mathcal{R}^{L}$, which have been optimized to distribute as tight as
possible, prepared for further processes, such as dimensionality reduction,
feature learning, and so on.

###  4.4. Discovery Algorithm in Latent Space

Figure 9. Example of Causal Discovery in feature space. Nodes indicate the
variable representations or latent factors.

Here we adopt heuristics to realize causal discovery in latent space among
Variable representations, using KLD (KL-Divergence) as the metric to evaluate
the similarity between two feature vectors in the common latent space
$\mathcal{R}^{L}$, as shown in Algorithm 1. The conventional MSE evaluation
can also be considered. However, it has been realized that the MSE-based
optimization subject to acyclicity constraint has some inherent restrictions
and may mislead the model to be dominated by data variance instead of
causality (reisach2021beware, ; kaiser2021unsuitability, ).

An illustrative example is given in Fig 9. In the four consecutive steps, edge
$e_{1}$ and edge $e_{3}$ are selected one after another, since the selected
$e_{1}$ make node B reachable, and B is the start node of $e_{3}$. At step
$3$, the learned causal effect of $e_{2}$, from $A$ to $C$, is deleted from
the candidate edges and need to be recalculated, because edge $e_{3}$ has
reached $C$ and changed the causal conditions on $C$.

Result: ordered edges set $\mathbf{E}=\\{e_{1},\ldots,e_{n}\\}$

$\mathbf{E}=\\{\\}$ ;

$N_{R}=\\{n_{0}\mid n_{0}\in N,Parent(n_{0})=\varnothing\\}$ ;

while  _$N_{R}\subset N$_ do

$\Delta=\\{\\}$ ;

for  _ $n\in N$ _ do

for  _$p\in Parent(n)$_ do

if  _$n\notin N_{R}$ and $p\in N_{R}$_ then

$e=(p,n)$;

$\beta=\\{\\}$;

for  _$r\in N_{R}$_ do

if  _$r\in Parent(n)$ and $r\neq p$_ then

$\beta=\beta\cup r$

end if

end for

$\delta_{e}=K(\beta\cup p,n)-K(\beta,n)$;

$\Delta=\Delta\cup\delta_{e}$;

end if

end for

end for

$\sigma=argmin_{e}(\delta_{e}\mid\delta_{e}\in\Delta)$;

$\mathbf{E}=\mathbf{E}\cup\sigma$;

$N_{R}=N_{R}\cup n_{\sigma}$;

end while

Algorithm 1 Causal Discovery with KLD Metric

$G=(N,E)$ | graph $G$ consists of nodes set $N$ and edges set $E$
---|---
$N_{R}$ | the set of reachable nodes
$\mathbf{E}$ | edges in order of being discovered
$K(\beta,n)$ | KLD metric between the causes set $\beta$ and effect node $n$
$\Delta=\\{\delta_{e}\\}$ | the set of all KLD Gain $\delta_{e}$ for each candidate edge $e$
$n$,$p$,$r$ | notation of node
$e$,$\sigma$ | notation of edge

##  5\. Feasibility Experiments

The feasibility experiments aim to confirm that: 1) the proposed autoencoder
with _Encrypt_ and _Decrypt_ can effectively extract higher-dimensional causal
variable representations. 2) the causal effects can be successfully
represented in latent space and used to complete graphical causal
reconstructions. 3) the heuristic causal discovery can effectively work in
latent space.

Our experiments are performed on a professionally calibrated synthetic
hydrology dataset generated by SWAT (Soil & Water Assessment Tool), a commonly
used hydrology data simulation system based on physical modules. These
physical models are jointly connected to form a causal graph. The question
from the hydrology domain is that the models learned from nearby watersheds
can hardly share their learned knowledge because existing physical methods
lack generalizability. The next experiment on this data will be toward a
practically effective model-generalization method.

However, this problem belongs to correlation but not necessarily to causality.
Due to some empirical restrictions, in this work, we do not have EHR data to
verify how well it works on CRB adjustment, but that would be our immediate
next move.

###  5.1. Hydrology Dataset

The hydrology system entails water fluxes and states across the earth’s space,
such as snowpack melting, evapotranspiration, and soil moisture. The
underlying hydrologic processes are highly complex, exhibit nonlinearities,
and contain a lot of unobservable dependencies. In recent years machine
learning on causal inferences naturally attracted hydrologists’ attention
(debate_Goodwill2020, ); Also, DL-based methods are widely used to extract
representations from time series efficiently. As a typical application, RNN
models contribute state-of-the-art techniques for streamflow prediction
(hess-18-Kratzert, ; Kratzert2019, ; 2022Li_Soil, ).

Table 1. Nodes in Fig 10 Hydrology Causal Graph

ID | Variable Name | Interpretation
---|---|---
A | Conditions Part I |  | Wind Speed $+$ Humidity $+$ Temperature
---
B | Conditions Part II |  | Temperature + Solar Radiation
---
\+ Precipitation
C | Evapotranspiration | Evaporation and transpiration
D | Snowpack | The winter frozen water in the ice form
E | Soil Water | Soil moisture in vadose zone
F | Aquifer | Groundwater storage
G | Surface Runoff | Flowing water over the land surface
H | Lateral | Vadose zone flow
I | Baseflow | groundwater discharge
J | Streamflow | Surface runoff + lateral flow + baseflow

Figure 10. The Underlying Causal Graph from domain knowledge with different
causality strengths.

_Streamflow_ is a critical state variable in hydrologic processes, where the
hydrologists attempted to explain the causal relationship between observable
environmental variables and the subsequential streamflow records, but have not
achieved a satisfying performance. Our experiments are based on such a
prediction task, using synthetic data generated by SWAT system, which
simulates the Root River Headwater watershed in Southeast Minnesota. From this
dataset, we adopted the time series that lasts 60 virtual years, in a daily
update frequency, to be the unsupervised training samples.

Fig 10 displays the hypothesized graphical causality assumed in SWAT.
Explanations about the nodes in graph are listed in Table 1. Hydrologic
processes present different levels of contributions. Surface runoff generation
processes (1st tier causality) cause streamflow peak more easily than lateral
flow generation processes (2nd tier causality), which are more important than
baseflow dynamics (3rd tier causality). For convenience, they are represented
by different colors.

###  5.2. Experimental Results of Causal Variable Representation Extraction

The proposed causal Variable representation extraction method is performed
sequentially on each node of the causal graph, from $A$ to $J$. The challenges
of these tasks come from the low dimensionality of variables. As shown in
Table 2, the number of attribute columns maximally reaches 5, and the
predictive target node $J$, i.e., Streamflow, is just a one-column variable.

We firstly expand the attributes of each node to a uniform larger count of 12
by randomly repeating their columns, then augmented with 12-dimensional dummy
variables of months to form a 24-dimensional input layer. In the 24 double-
wise feature amplifications, every column has been extracted 23 times,
augmented with the original vector. This way, we have a 576-dimensional output
through _Enrypt_ process. In our experimental context, the feature space for
representations is set to be 16-dimensional.

Table 2 shows each node’s statistics characteristics and the RMSE Evaluation
between the Output and Input of the autoencoder, where less RMSE indicates
lower information loss; in other words, the proposed autoencoder architecture
is more effectively extracted the causal variable’s representation. We provide
two RMSE evaluations for the nodes in the scaled (i.e., normalized) and
original values, respectively. Besides, the characteristic columns, including
Mean, Std, Min, and Max, are calculated at the scaled values.

It is worth mentioning that the hydrology dataset contains a large number of
meaningful zero values. For example, the variable of node $D$, named Snowpack,
only has non-zero values in winter, whose water body will present as Soil
Water in the other seasons. However, the zeros do have meanings since they
represent the variables vanishing. Therefore, we perform double
reconstructions simultaneously in autoencoder: one for the real value and the
other for the non-zero indicator (accessed by BCE), named Mask. The RMSE
performances in Table 2 are obtained by multiplying these two results. The
percentage of non-zero values is also provided for each node. These very low
RMSE values indicate the success of the reconstruction processes, which are in
the range of $[0.01,0.09]$, except node $F$, the Aquifer variable. It has been
known that the aquifer system modeling is still premature in the present
hydrology area. So this is reasonable to infer that in the synthetic data,
Aquifer is closer to random noise than other variables.

Table 2. Performance of Autoencoder (Causal Variables Representation).

Nodes | # Column | Mean | Std | Min | Max | Non-Zero Rate% | RMSE On Scaled | RMSE On Original | BCE On Mask
---|---|---|---|---|---|---|---|---|---
A | 5 | 1.8513 | 1.5496 | -3.3557 | 7.6809 | 87.54 | 0.093 | 0.871 | 0.095
B | 4 | 0.7687 | 1.1353 | -3.3557 | 5.9710 | 64.52 | 0.076 | 0.678 | 1.132
C | 2 | 1.0342 | 1.0025 | 0.0 | 6.2145 | 94.42 | 0.037 | 0.089 | 0.428
D | 3 | 0.0458 | 0.2005 | 0.0 | 5.2434 | 11.40 | 0.015 | 0.679 | 0.445
E | 2 | 3.1449 | 1.0000 | 0.0285 | 5.0916 | 100 | 0.058 | 3.343 | 0.643
F | 4 | 0.3922 | 0.8962 | 0.0 | 8.6122 | 59.08 | 0.326 | 7.178 | 2.045
G | 4 | 0.7180 | 1.1064 | 0.0 | 8.2551 | 47.87 | 0.045 | 0.81 | 1.327
H | 4 | 0.7344 | 1.0193 | 0.0 | 7.6350 | 49.93 | 0.045 | 0.009 | 1.345
I | 3 | 0.1432 | 0.6137 | 0.0 | 8.3880 | 21.66 | 0.035 | 0.009 | 1.672
J | 1 | 0.0410 | 0.2000 | 0.0 | 7.8903 | 21.75 | 0.007 | 0.098 | 1.088

Table 3. Brief Results of Heuristic Discovery: the edges in detected order and their measured causal strengths. Edge | A$\Rightarrow$C | B$\Rightarrow$D | C$\Rightarrow$D | C$\Rightarrow$G | D$\Rightarrow$G | G$\Rightarrow$J | D$\Rightarrow$H | H$\Rightarrow$J | B$\Rightarrow$E | E$\Rightarrow$G | E$\Rightarrow$H | C$\Rightarrow$E | E$\Rightarrow$F | F$\Rightarrow$I | I$\Rightarrow$J | D$\Rightarrow$I
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
KLD | 7.63 | 8.51 | 10.14 | 11.60 | 27.87 | 5.29 | 25.19 | 15.93 | 37.07 | 39.13 | 39.88 | 46.58 | 53.68 | 45.64 | 17.41 | 75.57
Gain | 7.63 | 8.51 | 1.135 | 11.60 | 2.454 | 5.29 | 25.19 | 0.209 | 37.07 | -5.91 | -3.29 | 2.677 | 53.68 | 45.64 | 0.028 | 3.384

###  5.3. Experimental Results of Causal Effect Representation and
Reconstruction

The causal effect representation learning in latent space is evaluated on each
pairwise causation. Then, by stacking them, we successfully reconstruct the
graphical causal structure shown in Fig 10. For evaluating the performances of
these stacking processes, say whether they can enhance the causal effect of
the _result_ node, we make such performance comparisons from the view of every
single node, to be as the effected _result_ node in a causal relation, as
listed in Table 4. For example, node $G$ has three possible causes $C$, $D$,
and $E$. Thus we can build up four models in total: the effect estimation of
each pairwise causation without stacking, including $C\Rightarrow G$,
$D\Rightarrow G$, $E\Rightarrow G$, and the effect estimation as in the fully
reconstructed graphical models, which is $CDE\Rightarrow G$. For convenience,
we call them the _pairwise-effect_ and _stacking-effect_ , respectively.
Additionally, in Table 4 we also display the initial performance of each node,
i.e., its autoencoder performance before learning any causal effects, to be
the baseline comparison, named as “Variable Representation (Before Effect
Learning)” in the table.

Figure 11.  Causal Graphical Structure Reconstruction Performance Comparision
(data simulations)

There are three different optimization tasks in the causal effect
representation learning process (refer to Section 4.3), among which both the
_second_ and _third_ one are about the _result_ node. Therefore, in Table 4 we
give out the performances of these two roles for each node: the role in the
_third_ optimization is shown as the columns named “Variable Representation
(After Effect Learning)”; and the role of the _second_ one is in the columns
named “Causal Effect Representation”, which represent the effect learning
performances of RNN models, along with the KLD metrics to present learned
causal strengthed.

By observations, the KLD metric differences in causal strength are easy to
see. Node $J$ has the smallest KLD values, which means that node $G$ (Surface
Runoff), $H$ (Lateral), and $I$ (Baseflow) significantly cause $J$
(Streamflow). While the weak causal relations present as high values of KLD.
For example, node $I$ can hardly be predictable given $D$ and $F$. For
_result_ nodes, $D$, $E$, $J$, the _stacking-effect_ causal strength is in the
intermediate range of _pairwise-effect_ strength, which implies that the
associative relationships among their causes nodes are ambiguous. For $G$ and
$H$, the _stacking-effect_ causal strength KLD is lower than all the
_pairwise-effect_ ones. It is reasonable to infer that the _stacking-effect_
model has captured additional causal information from the associative
relationships among their causes nodes. Besides, the KLD metric also tells
which cause node contributes the most to this causal effect. For example,
$C\Rightarrow G$ strength is closer to $CDE\Rightarrow G$ than the other
causes nodes, which indicates $C$ to be the strongest source in this causal
effect.

To better present the Causal Graph Structure Reconstruction performance, we
use these reconstructed models as a data simulator and show their simulation
results as in Fig 11, which displays three different nodes, $J$, $G$, and $I$
in a same synthetic year. For each node, the initial variable representation
performance (i.e., the autoencoder performance before effects learning) is
plotted as the blue line, and all the causal effects we have learned are
displayed in different colors. Besides RMSE, here we also adopted a hydrology-
using metric NSE (Nash–Sutcliffe model efficiency coefficient) to assess
predictive accuracy, where NSE $=1$ defines the best prediction. The three
initial variable representations (blue line) in Fig 11 almost reach the best
performance, as they overlapped with the ground truth observations (black
dots). The red line represents the simulations from _stacking-effect_ , whose
performances are highly consistent with our analysis above: Node $J$ has the
best prediction; node $I$ can hardly be predictable; and node G is predicted
best by $CDE$ causes combination, where $C$ provides the strongest causality.

###  5.4. Experimental Results of Causal Discovery

The discovery heuristic process is based on the KLD causal strength. 3 shows
the edges’ discovered order and the corresponding KLD and KLD-Gain. Here the
table cells are colored according to different tiers of causality, i.e., the
ground truth of causal strengths, as displayed in Fig 10. Since the discovery
algorithm successfully distinguished different causality strength levels, the
proposed KLD metric is confirmed effective.

Due to limited space, the result in Table 3 is a brief version, while the
complete one is provided in Appendix A Table 1, which includes all results in
each detecting round and the corresponding decision-making. We also performed
the conventional FGES (Fast Greedy Equivalence Search) causal discovery in
$10$ cross-validation, which is worse than our method, given in Appendix A
Table 2.

Table 4. The latent-space causal effect representations learning performances.

|  Variable Representation (Before Effect Learning) |  |  Variable Representation (After Effect Learning) |  Causal Effect Representation
---|---|---|---|---
|  RMSE | BCE |  |  RMSE | BCE |  RMSE | BCE | KLD
Effect Node |  | On Scaled
---
Values
| On Original
---
Values
Mask |  Causes Node(s) |  | On Scaled
---
Values
| On Original
---
Values
Mask |  | On Scaled
---
Values
| On Original
---
Values
Mask |  | In Feature
---
Space
C | 0.037 | 0.089 | 0.428 | A | 0.0295 | 0.0616 | 0.4278 | 0.1747 | 0.3334 | 0.4278 | 7.6353
|  |  |  | BC | 0.0350 | 1.0179 | 0.1355 | 0.0509 | 1.7059 | 0.1285 | 9.6502
|  |  |  | B | 0.0341 | 1.0361 | 0.1693 | 0.0516 | 1.7737 | 0.1925 | 8.5147
D | 0.015 | 0.679 | 0.445 | C | 0.0331 | 0.9818 | 0.3404 | 0.0512 | 1.7265 | 0.3667 | 10.149
|  |  |  | BC | 0.4612 | 26.605 | 0.6427 | 0.7827 | 45.149 | 0.6427 | 39.750
|  |  |  | B | 0.6428 | 37.076 | 0.6427 | 0.8209 | 47.353 | 0.6427 | 37.072
E | 0.058 | 3.343 | 0.643 | C | 0.5212 | 30.065 | 1.2854 | 0.7939 | 45.791 | 1.2854 | 46.587
F | 0.326 | 7.178 | 2.045 | E | 0.4334 | 8.3807 | 3.0895 | 0.4509 | 5.9553 | 3.0895 | 53.680
|  |  |  | CDE | 0.0538 | 0.9598 | 0.0878 | 0.1719 | 3.5736 | 0.1340 | 8.1360
|  |  |  | C | 0.1057 | 1.4219 | 0.1078 | 0.2996 | 4.6278 | 0.1362 | 11.601
|  |  |  | D | 0.1773 | 3.6083 | 0.1842 | 0.4112 | 8.0841 | 0.2228 | 27.879
G | 0.045 | 0.81 | 1.327 | E | 0.1949 | 4.7124 | 0.1482 | 0.5564 | 10.852 | 0.1877 | 39.133
|  |  |  | DE | 0.0889 | 0.0099 | 2.5980 | 0.3564 | 0.0096 | 2.5980 | 21.905
|  |  |  | D | 0.0878 | 0.0104 | 0.0911 | 0.4301 | 0.0095 | 0.0911 | 25.198
H | 0.045 | 0.009 | 1.345 | E | 0.1162 | 0.0105 | 0.1482 | 0.5168 | 0.0097 | 3.8514 | 39.886
|  |  |  | DF | 0.0600 | 0.0103 | 3.4493 | 0.1158 | 0.0099 | 3.4493 | 49.033
|  |  |  | D | 0.1212 | 0.0108 | 3.0048 | 0.2073 | 0.0108 | 3.0048 | 75.577
I | 0.035 | 0.009 | 1.672 | F | 0.0540 | 0.0102 | 3.4493 | 0.0948 | 0.0098 | 3.4493 | 45.648
|  |  |  | GHI | 0.0052 | 0.0742 | 0.2593 | 0.0090 | 0.1269 | 0.2937 | 5.5300
|  |  |  | G | 0.0077 | 0.1085 | 0.4009 | 0.0099 | 0.1390 | 0.4375 | 5.2924
|  |  |  | H | 0.0159 | 0.2239 | 0.4584 | 0.0393 | 0.5520 | 0.4938 | 15.930
J | 0.007 | 0.098 | 1.088 | I | 0.0308 | 0.4328 | 0.3818 | 0.0397 | 0.5564 | 0.3954 | 17.410

##  6\. Conclusion

In this paper, we initially discussed the concept of Multi-Dimensionality for
Acyclic Causal Graphs. It is the fundamental reason for inherent confounding
biases of causal learning, which we first named Causal Representation Bias
(CRB). In our opinion, this is also where the difficulties of applying the
Deep Learning models in causality stem from. Because CRB cannot be eliminated
from causal learning by reducing the modeling approximation error; On the
contrary, the highly efficient global optimization strategy with continuous
constraint may increase the influence of CRB.

In the discussed scenario, the 3-Dimensional space of DAG is defined by the
axes of the causal variables, the Absolute timeline where the causal learning
happens, and Relative timelines where the CRB possibly happens. We accordingly
proposed a novel new framework Causal Representation Learning (CRL), which is
to realize graphical models in latent space, such that 1) all causal relations
from the DAG can be transformed onto a common space to realize CRB adjustment
and de-confounded causal modeling; 2) the model interpretability can be
preserved in the latent space to realize Deep Learning causal models’
generalization.

We also proposed a novel deep network architecture (with the originally
designed layers Encrypt and Decrypt) to realize CRL and experimentally
confirmed its feasibility.

CRL provides a new perspective on causal learning and raises many valuable
questions as well as forward directions, like the meaning of causal graph
dimensionality, feature explanation in the latent learning space, and so on;
From the aspect of applications, like complex medical effects estimation, non-
model-based causal simulation, and so on. We expect a flourishing development
can happen in the following years.

## References

  * (1) Wood, Christopher J., and Robert W. Spekkens. _The lesson of causal discovery algorithms for quantum correlations: Causal explanations of Bell-inequality violations require fine-tuning._ New Journal of Physics 17.3 (2015): 033002.
  * (2) Vuković, Matej, and Stefan Thalmann. _Causal Discovery in Manufacturing: A Structured Literature Review._ Journal of Manufacturing and Materials Processing 6.1 (2022): 10.
  * (3) Ombadi, Mohammed, et al. _Evaluation of methods for causal discovery in hydrometeorological systems._ Water Resources Research 56.7 (2020): e2020WR027251.
  * (4) Zheng, Xun, et al. _Learning sparse nonparametric dags._ International Conference on Artificial Intelligence and Statistics. PMLR, 2020.
  * (5) Lachapelle, Sébastien, et al. _Gradient-based neural dag learning._ arXiv preprint arXiv:1906.02226 (2019).
  * (6) Luo, Yunan, Jian Peng, and Jianzhu Ma. _When causal inference meets deep learning._ ” Nature Machine Intelligence 2.8 (2020): 426-427.
  * (7) Sobel, Michael E. _An introduction to causal inference._ Sociological Methods & Research 24.3 (1996): 353-379.
  * (8) Glymour, Clark, Kun Zhang, and Peter Spirtes. _Review of causal discovery methods based on graphical models._ Frontiers in genetics 10 (2019): 524.
  * (9) Li, Jia, et al. _Teaching deep learning causal effects improves predictive performance._ arXiv preprint arXiv:2011.05466 (2020).
  * (10) Schölkopf, Bernhard, et al. _Toward causal representation learning._ Proceedings of the IEEE 109.5 (2021): 612-634.
  * (11) Li, Yunzhu, et al. _Causal discovery in physical systems from videos._ Advances in Neural Information Processing Systems 33 (2020): 9180-9192.
  * (12) Zhong, Guoqiang, et al. _An overview on data representation learning: From traditional feature learning to recent deep learning._ The Journal of Finance and Data Science 2.4 (2016): 265-278.
  * (13) Bengio, Yoshua, Aaron C. Courville, and Pascal Vincent. _Unsupervised feature learning and deep learning: A review and new perspectives._ CoRR, abs/1206.5538 1.2665 (2012): 2012.
  * (14) Goodwell, Allison E., et al. _Debates—Does information theory provide a new paradigm for Earth science? Causality, interaction, and feedback._ Water Resources Research 56.2 (2020): e2019WR024940.
  * (15) Kratzert, Frederik, et al. _Rainfall–runoff modeling using long short-term memory (LSTM) networks._ Hydrology and Earth System Sciences 22.11 (2018): 6005-6022.
  * (16) Kratzert, Frederik, et al. _Towards learning universal, regional, and local hydrological behaviors via machine learning applied to large-sample datasets._ Hydrology and Earth System Sciences 23.12 (2019): 5089-5110.
  * (17) Li, Xiaojie, Jian Tang, and Chengxiang Yin. _Sequence-to-Sequence Learning for Prediction of Soil Temperature and Moisture._ IEEE Geoscience and Remote Sensing Letters 19 (2022): 1-5.
  * (18) Reisach, Alexander G., Christof Seiler, and Sebastian Weichwald. _Beware of the simulated dag! varsortability in additive noise models._ arXiv preprint arXiv:2102.13647 (2021).
  * (19) Kaiser, Marcus, and Maksim Sipos. _Unsuitability of NOTEARS for causal graph discovery._ arXiv preprint arXiv:2104.05441 (2021).
  * (20) Zheng, Xun, et al. _Dags with no tears: Continuous optimization for structure learning._ Advances in Neural Information Processing Systems 31 (2018).
  * (21) Dinh, Laurent, Jascha Sohl-Dickstein, and Samy Bengio. _Density estimation using real nvp._ arXiv preprint arXiv:1605.08803 (2016).

##  Appendix A Complete Experimental Results of Causal Discovery

See pages - of mainz_append_pdf.pdf


# Mass Superinflation from Enhanced Cauchy Horizon Singularity
in a Reissner-Nordström Black Hole

Nihar Ranjan<EMAIL_ADDRESS>and Malay K.
<EMAIL_ADDRESS>
Department of Physics, Indian Institute of Technology Guwahati
Guwahati 781 039, India

(January 20, 2024)

###### Abstract

Ever since Penrose and Simpson contradicted Novikov’s prediction that an
infalling passenger would emerge into an asymptotically flat universe, there
have been a continued interest in predicting the nature of singularity at the
Cauchy horizon of a Reissner-Nordström black hole. This prediction was first
confirmed by Poisson and Israel using cross-stream of massless particles,
suggesting the phenomenon of mass inflation. Ori however obtained a weaker
singularity using a null shell of radiation. In this work, we consider a
massive scalar field coupled to the Reissner-Nordström geometry and analyze
the nature of singularity at the Cauchy horizon. To study the asymptotic
behavior of the mass function and the scalar field near the Cauchy horizon, we
perturbatively solve the coupled dynamical equations emplyoing the Adomian
decomposition method. Our analysis shows that the mass function exhibits a
very rapid and unbounded double-exponential growth, called herein mass
superinflation, which is enormously stronger than previously obtained
singularities. The scalar field is also found to undergo a very strong
blueshift near the Cauchy horizon.

##  1 Introduction

The result of any general gravitational collapse is a trapped region, which
indicates the presence of both an outer and an inner apparent horizons.
Classically the outer apparent horizon shifts or moves outward (as more mass
falls into the black hole), while the inner one moves inwards, either settling
to a final constant radius or collapsing all the way to the origin and
creating a spacelike singularity. If the inner horizon tends to a final
position of constant radius then it leads to the well-known mass inflation
instability phenomenon, where any perturbation in the matter content of the
geometry leads to the increase of curvature at the Cauchy horizon (CH). The
most common geometries with an inner horizon are the Kerr and Reissner-
Nordström solutions. In this work we study the behaviour at the CH in the
Reissner-Nordström (RN) black hole.

A Reissner-Nordström spacetime is formed when an electrically charged
spherically symmetric mass undergoes sufficient gravitational collapse. The
static geometry of the spacetime is obtained as a unique solution of the
Einstein-Maxwell equation [1],[2] and the spacetime metric is described by

$ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}d\Omega^{2}~{}.$ (1)

where $f(r)=(1-\frac{2m_{0}}{r}+\frac{Q^{2}}{r^{2}})$ is the redshift
function, with $m_{0}$ the ADM mass and $Q$ the electric charge of the black
hole (BH) and $d\Omega^{2}=d\theta^{2}+\sin^{2}\theta d\phi^{2}$. In this
geometry, there are two horizons, given by $f(r)=0$,

$r_{\pm}=m_{0}\pm\sqrt{m_{0}^{2}-Q^{2}}~{}.$

The inner horizon $(r_{-})$ of the Reissner-Nordström black hole also serves
as the Cauchy Horizon, because the zero of the redshift function has a
negative slope at the inner horizon. Hence the inner horizon acts as a
boundary for predictability. The cosmic censorship conjecture implies
$m_{0}\geq Q$, with the equality (inequality) giving extremal (non-extremal)
black hole [3].

Figure 1: Penrose diagram of a Reissner-Nordström black hole showing the
trajectory of an infalling passenger who manages to skip the central
singularity and emerges into another universe.

The extended Penrose diagram for non-extremal Reissner-Nordström black hole
with an infalling passenger in this background is shown in Figure 1. The
redshift function, $f(r)$, changes sign in the region $r_{-}<r<r_{+}$, which
implies that the roles of $t$ and $r$ coordinates are interchanged.

The presence of a timelike singularity gives rise to a strange phenomenon.
Novikov [4] first argued that an object falling into an RN black hole can
emerge into an asymptotically flat universe without hitting the central
singularity. During gravitational collapse with a given initial condition, the
internal dust evolves in a completely regular manner, creating a tunnel to
another universe, as shown in Figure 1.

However, this deterministic behavior is avoided by the instability of the
inner horizon, serving as a barrier to predictability. Simpson and Penrose [5,
6], shortly after the formulation of the singularity theorem, considered an
infalling test electromagnetic field in the RN spacetime. They confirmed that
the (CH) becomes unstable as the test field develops a singularity. Thus an
infalling passenger will not survive upon encountering the CH.

Penrose’s prediction of CH instability can be understood by considering an
infalling test field [7] such as a wave with uniformly spaced crests, which
will become densely packed and experience gravitational and Doppler
blueshifts. This amplification continues indefinitely at the CH signifying its
instability to any external perturbation [8]. Owing to the backreaction of
this divergence, the regular CH transforms into a curvature singularity. After
Penrose, several confirmations were made regarding this divergence up to
linear perturbation, but no analysis was carried out beyond linear order until
Poisson and Israel.

###  1.1  Classical Mass Inflation

Poisson and Israel [9], first examined the behavior of the mass function
accounting for the back reaction on the CH in the RN geometry due to
infinitely blueshifted infalling spherically symmetric stream of massless
particles. This resulted in a charged Vaidya geometry and, for a suitable
choice of coordinates, the singularity was weak. However, upon including the
backscattered outgoing stream, the mass function behaved as $m(U,V)\sim
m_{0}\epsilon^{2}\delta(U)e^{k_{0}v}(v/m_{0})^{-12}$. Thus within a time scale
of $\Delta v\sim k_{0}^{-1}\sim m_{0}$, the curvature $m/r^{3}$ must reach
Planckian values near the CH (at constant $U$), where $v$ is the retarted
time, $\epsilon$ is the quadrupole moment, $\delta(U)$ is the fraction of mass
radiated outward between the moment of entry and time $U$. The outgoing
radiation is the main reason behind this sudden growth in strength of the
singularity, as was first found by Poisson [10].

In this process the mass function diverges exponentially (coined as mass
inflation) and the Kretschmann scalar
$R_{\mu\nu\alpha\beta}R^{\mu\nu\alpha\beta}$ approaches infinity at the CH.
This is a stronger singularity than the charged Vaidya solution. The reason
behind this strong singularity is the failure of the cross streams to move
faster than light with respect to each other, which produces a large amount of
gravitational energy that reveals itself as an infinite amount of local mass
near the CH [11].

Ori [7] improved the analysis by patching two Vaidya solutions separated by an
infalling null shell of radiation, yielding $m(v)\approx\Delta
m(v)\propto|v|^{-1}(-\ln|v|)^{-p}$, where $p\geq 12$ from Price’s law[12].
Thus $m(v)$ and $R_{\mu\nu\alpha\beta}R^{\mu\nu\alpha\beta}\propto m^{2}$
diverge near the CH.

Ori showed that the growth rate of curvature and tidal forces is proportional
to $\tau^{-2}|\ln|\tau||^{-p}$, $\tau$ being the proper time and particle
separations constituting an infalling passenger remain finite as
$\tau\rightarrow 0$, so that the spacetime can be extended beyond the CH
despite the (weak) mass inflation singularity. Singularity strength criteria
by Tipler [13], Ellis and Schmidt [14], and others [15, 16] confirm the same
conclusion.

Consequently attempts have been made in semiclassical frameworks to address
this problem.

###  1.2 Semiclassical Mass Inflation

Employing the semiclassical treatment of Davies and Fulling [17] for the
stress energy tensor, Hiscock [18] first showed that the stress energy tensor
diverges at the Cauchy horizon. Moreover, Balbinot and Poisson [19] considered
the semiclassical effect [20] on the mass inflation including quantum
corrections to the stress energy tensor. Although they did not give a definite
conclusion, they speculated that the mass inflation singularity may become
stronger or regular depending on the sign of the quantum correction.

Hwang and Yeom [21] examined the interior structure during the formation and
evaporation of RN black holes in the semiclassical formulation. Employing
S-wave approximation they obtained a wormhole-like internal structure when
$Q/M$ is small. However, in the final stages, mass inflation generates a large
curvature in the interior.

Recently Barceló [22, 23] analyzed the backreaction of a quantum field on the
RN geometry. They incorporated, in the Polykov approximation, the renormalized
stress energy tensor of a massless quantized scalar field and found that the
classical tendency of the inner horizon to move inward may be overcome by the
semiclassical tendency to move outward.

Carballo-Rubio et al. [24] and Franzin [25] took regularized black holes and
argued that they have vanishing inner horizon surface gravity. However, Maken
[26] proved that the CH instability still continues to exist.

###  1.3 de Sitter Models

Motivated by the accelerated expansion of the Universe, Brady and Poisson [27]
modeled the RN BH in de Sitter space, perturbed by an infalling stream of
particles. Consequently the CH instability occurs when its surface gravity
exceeds that of the cosmological horizon. However, for a generic perturbation,
the CH remains stable, requiring $|Q|>M$ [27, 28].

Brady et al. [29] showed for a perturbed RN BH that for specific region in the
parameter space the mass function remains bounded while the curvature scalar
diverges.

Cai and Su [30] employed the stability conjecture of Helliwell and Konkowski
[31, 32], confirming the prediction of Brady et al. However, Marković and
Poisson [33] claimed that the CH must be unstable when treated quantum
mechanically.

###  1.4 Present Work

We see from the above discussions that the problem of mass inflation has been
attacked from various angles in both classical and semiclassical frameworks.
All these studies lead to one common feature that the singularity of the CH
continues to exist both classically and semiclassically. One important feature
which remains still unclear is that whether the CH instability is strong
enough so that the possibility of extending the spacetime beyond the CH does
not exist. Another common feature that can be noticed is that those
calculations were based on ingoing and outgoing null shells of massless
radiation. A more physical approach would be to consider a simple but
realistic model describing the scenario by a massive scalar field as
astrophysical black holes are normally surrounded by massive particles in
addition to massless radiation.

In this paper, we therefore consider a massive, chargeless classical scalar
field in spherically symmetric Reissner-Nordström geometry and study the
behaviour of the singularity at the Cauchy horizon. We find that the
singularity is enhanced double exponentially as the CH is approached
asymptotically, where the scalar field undergoes an extraordinary blueshift.
These asymptotic behaviours are much stronger than those in previous studies
and they suggest that the inner horizon serves as an impenetrable physical
boundary of the spacetime.

The rest of the paper is organised as follows. In Section 2, we introduce our
model of massive scalar field in the Reissner-Nordström geometry and obtain
the Einstein-Maxwell field equations coupled with dynamics of the massive
scalar field. In Section 3, we develop a perturbation theory to solve the
coupled set of differential equations which gives the solutions in terms of an
integro-differential equation. In Section 4, we give a brief account of the
Adomian decomposition method (ADM) to be utilised in Section 5, where we
employ this method to solve the equation obtained in Section 3 in the
asymptotic limit. Finally in Section 6, we give a discussion and conclusion of
our results.

##  2 The Model

As discussed above, we study the behaviour of the CH singularity by perturbing
the black hole with a massive, chargeless, scalar field $\Phi$. Thereby, we
construct the dynamics of the scalar field coupled with the dynamics of the
mass function. In order to describe mass inflation, we shall solve the coupled
dynamical equations in the near horizon and long time limits.

We begin our calculation by modeling the phenomenon with the combined action

$S=S_{g}+S_{\Phi}~{},$ (2)

with

$S_{g}=\int
d^{4}x\sqrt{-g}\left[\frac{1}{2}M_{P}^{2}R-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}\right]~{},$

and

$S_{\Phi}=\int
d^{4}x\sqrt{-g}\left[-\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\Phi\nabla_{\nu}\Phi-\frac{1}{2}\mu^{2}\Phi^{2}\right]$

where $R$ is the Ricci scalar, $\mu$ is the mass of the scalar field $\Phi$
and $F_{\mu\nu}$ is the electromagnetic tensor. The spacetime metric
$g_{\mu\nu}$ is spherically symmetric in the Reissner-Nordstöm geometry and it
is given by the line element

$ds^{2}=f[-dt^{2}+dr^{*2}]+r^{2}d\Omega^{2}~{},$ (3)

where the tortoise coordinate $r^{*}$ is defined by $dr^{*}=dr/f(r)$. For
further convenience, we rewrite the line element in $(v,r)$ coordinateS as

$ds^{2}=-fdv^{2}+2dvdr+r^{2}d\Omega^{2}~{}.$ (4)

Now, because of the perturbation by the massive scalar field, the mass $m_{0}$
will be modified to $m_{0}+m(v)$, so that

$f=1-\frac{2[m_{0}+m(v)]}{r}+\frac{Q^{2}}{r^{2}}~{}.$

Henceforth, we shall denote the inner horizon ($r_{-}$) by $r_{0}$.

Equipped with the above formulation, we obtain from the action (2), the
Einstein field equation,

$\displaystyle\begin{split}R_{\mu\nu}-&\frac{1}{2}g_{\mu\nu}R=\frac{1}{M_{p}^{2}}\left[F_{\nu\beta}F_{\mu}^{\beta}-\frac{1}{4}g_{\mu\nu}F_{\alpha\beta}F^{\alpha\beta}\right.\\\
&\left.+\nabla_{\mu}\Phi\nabla_{\nu}\Phi-
g_{\mu\nu}\left(\frac{1}{2}\nabla^{\rho}\Phi\nabla_{\rho}\Phi+\frac{1}{2}\mu^{2}\Phi\right)\right]~{}.\end{split}$
(5)

The $00$ component of the above equation yields the dynamics of the mass
function as

$\begin{split}\frac{dm}{dv}=4\pi
r^{2}&\left[\left(\frac{\partial\Phi}{\partial
v}\right)^{2}+f\left(\frac{\partial\Phi}{\partial
v}\right)\left(\frac{\partial\Phi}{\partial r}\right)\right.\\\
&+\left.\frac{1}{2}f^{2}\left(\frac{\partial\Phi}{\partial
r}\right)^{2}+\frac{1}{2}f\mu^{2}\Phi^{2}\right]~{}.\end{split}$ (6)

Moreover, the action (2) gives the Klein-Gordon equation for the massive
scalar field, $\left[\square-\mu^{2}\right]\Phi=0$, which can be written as

$\begin{split}2r^{2}\frac{\partial^{2}\Phi}{\partial v\partial
r}+2r\frac{\partial\Phi}{\partial v}&+2rf\frac{\partial\Phi}{\partial
r}+r^{2}\frac{\partial f}{\partial r}\frac{\partial\Phi}{\partial r}\\\
&+r^{2}f\frac{\partial^{2}\Phi}{\partial r^{2}}-\mu^{2}\Phi
r^{2}=0~{}.\end{split}$ (7)

Equation (7), giving the dynamics of the massive scalar field $\Phi$, is
coupled to equation (6) describing how the mass function $m(v)$ of the black
hole behaves.

In order to study the asymptotic behaviour of the mass function $m(v)$ near
the CH, we need to solve the two coupled, inhomogeneous, nonlinear, partial
differential equations, given by equations (6) and (7) in the near horizon and
long time limit.

##  3 Asymptotic Dynamics

We can simplify the coupled equations by assuming that the massive scalar
field, $\Phi(r,v)$, is spherically symmetric (S wave approximation) so that
its intensity falls off like $1/r^{2}$, and we write

$\Phi(r,v)=\frac{1}{r}\phi(r,v)$ (8)

Substituting equation (8) in equations (6) and (7), we obtain

$\begin{split}&\frac{dm}{dv}=4\pi\left[\left(\frac{\partial\phi}{\partial
v}\right)^{2}+\frac{1}{2}\mu^{2}f\phi^{2}\right.\\\ &\hskip
42.67912pt\left.+f\frac{\partial\phi}{\partial
v}\left(\frac{\partial\phi}{\partial
r}-\frac{\phi}{r}\right)+\frac{1}{2}f^{2}\left(\frac{\partial\phi}{\partial
r}-\frac{\phi}{r}\right)^{2}\right]\\\ \end{split}$ (9)

and

$\begin{split}&rf\frac{\partial^{2}\phi}{\partial
r^{2}}+2r\frac{\partial^{2}\phi}{\partial v\partial r}+r\frac{\partial
f}{\partial r}\frac{\partial\phi}{\partial r}-\phi\frac{\partial f}{\partial
r}-\mu^{2}\phi r=0~{}.\\\ \end{split}$ (10)

To explore the dynamics near the CH, we expand $f(r,v)$ and $\phi(r,v)$ about
the inner horizon, defined by $r_{0}=m_{0}-\sqrt{m_{0}^{2}-Q^{2}}$,

$\begin{split}f(r,v)&=f_{0}(v)+xf_{1}(v)+x^{2}f_{2}(v)+\dots\\\
\phi(r,v)&=\phi_{0}(v)+x\phi_{1}(v)+x^{2}\phi_{2}(v)+\dots\\\ \end{split}$
(11)

where $x=r-r_{0}$. In the near horizon approximation $x\xrightarrow[]{}0$, we
can neglect the contributions from $f_{2}$, as it is expected to be much
smaller compared to the first-order contribution from $f_{1}$. The important
point to note here is that we are only interested in obtaining the behaviour
of the black hole when the infalling ovserver is in between the two horizons,
where $f(r,v)$ is negative. Thus

$f_{0}(v)=\frac{2m(v)}{r_{0}}$

and

$f_{1}(v)=-\frac{2}{r_{0}^{2}}\left[m_{0}+m(v)-\frac{Q^{2}}{r_{0}}\right]~{}.$

Now substituting the series expansions from equations (11) into equations (9)
and (10) yeilds

$\begin{split}\frac{dm}{dv}&=4\pi\left[\left(\frac{d\phi_{0}}{dv}\right)^{2}+\frac{1}{2}f_{0}\mu^{2}\phi_{0}^{2}\right.\\\
&+f_{0}\frac{d\phi_{0}}{dv}\left(\phi_{1}-\frac{\phi_{0}}{r_{0}}\right)\left.+\frac{1}{2}f_{0}^{2}\left(\phi_{1}-\frac{\phi_{0}}{r_{0}}\right)^{2}\right]\end{split}$
(12)

and

$\frac{d\phi_{1}}{dv}+f_{0}\phi_{2}+\frac{f_{1}}{2}\left(\phi_{1}-\frac{\phi_{0}}{r_{0}}\right)-\frac{1}{2}\mu^{2}\phi_{0}=0~{},$
(13)

at zeroth order in $x$. And from equation (10)

$\begin{split}\frac{d\phi_{2}}{dv}+f_{1}\phi_{2}=\frac{1}{4}\left[\mu^{2}\phi_{1}+\frac{f_{1}}{r_{0}}\left(\phi_{1}-\frac{\phi_{0}}{r_{0}}\right)\right]\end{split}$
(14)

at first order in $x$, which can be re-expressed as

$\begin{split}4r_{0}&\frac{d}{dv}\left[\left(\phi_{2}+\frac{\phi_{0}}{4r_{0}^{2}}\right)\exp\left(\int
f_{1}dv\right)\right]\\\
&=e^{-\mu^{2}r_{0}v}\frac{d}{dv}\left[\phi_{1}\exp\left(\int(\mu^{2}r_{0}+f_{1})dv\right)\right]\\\
&-\exp\left(\int
f_{1}dv\right)\frac{d}{dv}\left(\phi_{1}-\frac{\phi_{0}}{r_{0}}\right)~{}.\end{split}$

In the Cauchy horizon limit, $v\xrightarrow[]{}\infty$, the first term on the
right hand side vanishes, and we obtain

$\begin{split}4r_{0}&\frac{d}{dv}\left[\left(\phi_{2}+\frac{\phi_{0}}{4r_{0}^{2}}\right)\exp\left(\int
f_{1}dv\right)\right]\\\ &=-\exp\left(\int
f_{1}dv\right)\frac{d}{dv}\left(\phi_{1}-\frac{\phi_{0}}{r_{0}}\right)~{}.\end{split}$
(15)

The nonlinearity of the original equations (9) and (10) introduces a closure
problem in the above perturbed equations. This closure problem can be
circumvented by balancing equation (15) with the approximations
$\phi_{2}\approx-\phi_{0}/(4r_{0}^{2})$ and $\phi_{1}\approx\phi_{0}/r_{0}$.
These approximations are expected to work well for the near extremal case
$m_{0}\sim Q$. Thus the series expansion for $\phi(r,v)$ takes the form

$\phi(r,v)\bigg{|}_{\text{near
CH}}=\phi_{0}\left[1+\frac{x}{r_{0}}-\frac{x^{2}}{4r_{0}^{2}}+\dots\right]~{}.$
(16)

With these approximations, equation (13) gives

$\begin{split}&\frac{d\phi_{0}}{dv}=\frac{1}{2}\Biggl{(}r_{0}\mu^{2}+\frac{m(v)}{r_{0}^{2}}\Biggl{)}\phi_{0}\\\
\end{split}$ (17)

in the CH limit. Solution of this equation is readily obtained as

$\phi_{0}(v)=A\phantom{0}e^{\mu^{2}r_{0}v/2}\exp\bigg{(}\frac{1}{2r_{0}^{2}}\int
mdv\bigg{)}~{},$ (18)

where $A$ is the integration constant. Substituting this solution in the mass
equation (12) leads to

$\begin{split}\frac{r_{0}^{4}}{A^{2}\pi}&\exp\left(-\int\frac{m}{r_{0}^{2}}dv\right)\frac{dm}{dv}\\\
&=e^{\mu^{2}r_{0}v}\big{(}m^{2}+6\mu^{2}r_{0}^{3}m+r_{0}^{6}\mu^{4}\big{)}\end{split}$
(19)

This last mass equation (19) is in the form of an integro-differential
equation. For the sake of convenience, we reexpress it as a second order
differential equation of the form

$\begin{split}\frac{d^{2}h}{dv^{2}}=e^{h+av}\biggl{[}\bigg{(}\frac{dh}{dv}\bigg{)}^{2}+6a\frac{dh}{dv}+a^{2}\biggl{]}\end{split}$
(20)

with $a=\mu^{2}r_{0}$, $A^{2}=r_{0}^{4}/\pi$, and

$h(v)=\frac{1}{r_{0}^{2}}\int m(v)\,dv,$

so that

$m=r_{0}^{2}\frac{dh}{dv}~{}.$

To solve the second order differential equation (20), we employ the standard
Adomian decomposition method (ADM) with the initial condition that
$h(v_{0})=\alpha>0$ and $h^{\prime}(v_{0})=\beta>0$, where $v_{0}$ is assumed
to be the moment of crossing the outer horizon.

##  4 Adomian Decomposition Method

In this section, we give a brief account of the Adomian decomposition method
(ADM). For more detailed discussions, we invite the reader to Refs. [34, 35,
36, 37, 38], and references therein.

The ADM is a perturbative method, useful for solving a variety of ordinary and
partial linear or nonlinear (homogeneous or inhomogeneous) differential
equations. The solution for the dependent variable,
$y(x_{1},x_{2},\dots,x_{n})$ is given in terms of an infinite series, which
converges if an analytic solution exists.

Consider the general inhomogeneous, nonlinear ordinary differential equation,

$L[y]+N[y]+R[y]=g(x)~{},$ (21)

where $L$ is the highest order linear differential operator, $N$ is the
nonlinear part of the equation and $R$ is the remaining linear part, and $g$
is the source. If $L=\frac{d^{k}}{dx^{k}},$ assuming that the inverse of $L$
exists, then $L^{-1}$ is defined as

$L^{-1}=\int_{x_{0}}^{x}dx\int_{x_{0}}^{x}dx\dots\int_{x_{0}}^{x}dx\phantom{000}(\text{$k$
times})~{},$

where $x_{0}$ is the initial value of $x$.

Applying $L^{-1}$ on equation (21) yeilds

$L^{-1}L[y]=L^{-1}g-L^{-1}N[y]-L^{-1}R[y]~{}.$ (22)

The ADM prescribes that the dependent variable $y$ can be written as an
infinite sum of the form

$y(x)=\sum_{n=0}^{\infty}y_{n}(x)~{},$ (23)

and the nonlinear part $N[y]$ can be decomposed in terms of a series of
polynomials, known as Adomian polynomials, such that

$\begin{split}N[y]=\sum_{n=0}^{\infty}A_{n}~{},\end{split}$

with

$A_{n}=\frac{1}{n!}\frac{d^{n}}{d\lambda^{n}}\left[N\left(\sum_{p=0}^{\infty}y_{p}\lambda^{p}\right)\right]\Bigg{|}_{\lambda=0}~{}.$
(24)

Putting these summations in equation (22) leads to

$\begin{split}L^{-1}L\sum_{n=0}^{\infty}y_{n}(x)=L^{-1}g-&L^{-1}\bigg{(}R\sum_{n=0}^{\infty}\big{[}y_{n}(x)\big{]}\bigg{)}\\\
&-L^{-1}\bigg{(}\sum_{n=0}^{\infty}A_{n}\bigg{)}~{}.\end{split}$ (25)

In our case, $L$ is of second order $L=\frac{d^{2}}{dx^{2}}$, and we get the
recurrence relations

$\begin{split}&y_{0}=y(x_{0})+y^{\prime}(x_{0})(x-x_{0})+L^{-1}g\\\
&y_{1}=-L^{-1}(y_{0}))-L^{-1}R(y_{0})-L^{-1}(A_{0})\\\ &\dots,\\\
&y_{n+1}=-L^{-1}(y_{n})-L^{-1}R(y_{n})-L^{-1}(A_{n})~{}.\end{split}$ (26)

In some cases, the nonlinear part may have a very complicated structure and
calculation of the Adomian polynomials $(A_{n})$ become very cumbersome. In
such cases, $A_{n}$ can be calculated using a Taylor series. Supposing the
nonlinear part is some function of $y$ denoted by $F(y)$, the Taylor expansion
of $F(y)$ around $y_{0}$ is given by

$\begin{split}F(y)=F(y_{0})&+F^{\prime}(y_{0})(y-y_{0})\\\
&+\frac{1}{2!}F^{\prime\prime}(y_{0})(y-y_{0})^{2}+\dots~{}.\end{split}$

Upon substituting the series expansion (23) for $y(x)$ in the above Taylor
expansion yields

$\begin{split}F(y)=F(y_{0})&+F^{\prime}(y_{0})(y_{1}+y_{2}+...)\\\
&+\frac{1}{2!}F^{\prime\prime}(y_{0})(y_{1}+y_{2}+...)^{2}+...\end{split}$

The required Adomian polynomials can be obtained from the above Taylor
expansion as

$\begin{split}&A_{0}=F(y_{0})\\\ &A_{1}=F^{\prime}(y_{0})y_{1}\\\
&A_{2}=F^{\prime}(y_{0})y_{2}+\frac{1}{2!}F^{\prime}(y_{0})y_{1}^{2}\\\
&A_{3}=F^{\prime}(y_{0})y_{3}+\frac{2}{2!}F^{\prime\prime}(y_{0})y_{1}y_{2}+\frac{1}{3!}F^{\prime\prime\prime}(y_{0})y_{1}^{3}\end{split}$
(27)

We shall utilise this methodology in the next section to solve our
differential equation (20).

##  5 Solution using ADM

It is easy to see from equations (20) and (21) that $R=0$ and $g=0$ in our
present case. The nonlinear part is given by the function

$F=e^{h+av}[(h^{\prime})^{2}+6ah^{\prime}+a^{2}]$

and

$L=\frac{d^{2}}{dv^{2}}~{}\text{giving}~{}L^{-1}=\int_{v_{0}}^{v}dv\int_{v_{0}}^{v}dv.$

As stated earlier, we have the initial conditions $h[v_{0}]=\alpha$ and
$h^{\prime}[v_{0}]=\beta=m_{0}/r_{0}^{2}$, and the mass function $m(v)$ is
related to $h(v)$ by

$m(v)=r_{0}^{2}\frac{dh}{dv}=r_{0}^{2}\frac{d}{dv}\left[\sum_{n=0}^{\infty}h_{n}(v)\right].$

Thus from equation (20), we have

$\begin{split}L^{-1}L[h]=L^{-1}N[h]~{},\end{split}$

giving

$h[v]=h[v_{0}]+(v-v_{0})h^{\prime}[v_{0}]+L^{-1}\left[\sum_{n=0}^{\infty}A_{n}\right]~{}.$

We therefore obtain $h_{0}(v)=\alpha+\beta(v-v_{0})$, and the remaining
$h_{n}$’s are given by the recurrence relation

$h_{n+1}(v)=L^{-1}[A_{n}]\phantom{0},\hskip 28.45274ptn\geq 0$

where the Adomian polynomials, $A_{n}$, can be calculated using equation (27).
Some of the first few $h_{n}(v)$ are listed below:

$\begin{split}h_{0}(v)&=\alpha+\beta(v-v_{0})\\\ h_{1}(v)&=e^{(\alpha-\beta
v_{0})}(a^{2}+6a\beta+\beta^{2})\left[e^{v(a+\beta)}\right.\\\
&\left.+e^{v_{0}(a+\beta)}\left(-1-(v-v_{0})(a+\beta)\right)\right]\\\
h_{2}(v)&=\frac{e^{2(\alpha-\beta
v_{0})}(a^{2}+6a\beta+\beta^{2})^{2}}{4\beta(a+\beta)^{3}}\left[e^{2v(a+\beta)}\right.\\\
&+e^{(v+v_{0})(a+\beta)}\\{4-4(v-v_{0})(a+\beta)\\}\\\
&\left.+e^{2v_{0}(a+\beta)}\\{-5-2(v-v_{0})(a+\beta)\\}\right]\\\
h_{3}(v)&=\frac{e^{3(\alpha-\beta
v_{0})}(a^{2}+6a\beta+\beta^{2})^{3}}{12\beta^{2}(a+\beta)^{4}}\Biggl{[}e^{3v(a+\beta)}\\\
&+e^{3v_{0}(a+\beta)}\\{-22-6(v-v_{0})(a+\beta)\\}\\\
&+e^{(2v+v_{0})(a+\beta)}\\{6-6(v-v_{0})(a+\beta)\\}\\\
&+3e^{(v+2v_{0})(a+\beta)}\Biggl{\\{}5+2(v-v_{0})(a+\beta)\\\
&\times\bigg{(}-3+(v-v_{0})(a+\beta)\bigg{)}\Biggl{\\}}\Biggl{]}\\\
\end{split}$ (28)

The corresponding mass function is readily obtained as

$\begin{split}\frac{m(v)}{r_{0}^{2}}&=\Bigg{[}h_{o}^{\prime}(v)+h_{1}^{\prime}(v)+h_{2}^{\prime}(v)+h_{3}^{\prime}(v)\dots\Bigg{]}\\\
&=\beta+\frac{1}{\sigma}\Bigg{[}\rho\sigma(v-v_{0})e^{(\delta+v)}\\\
&\phantom{0000}+\frac{\rho^{2}\sigma^{2}}{6}(v-v_{0})^{3}e^{2(\delta+v)}\\\
&\phantom{0000}+\frac{\rho^{3}\sigma^{3}}{30}(v-v_{0})^{5}e^{3(\delta+v)}+\dots\Bigg{]}~{},\end{split}$

which can be written as

$\frac{m}{r_{0}^{2}}=\beta+\frac{1}{\sigma}+\frac{1}{\sigma(v-v_{0})}\sum_{n=0}^{\infty}\frac{z^{n}}{n!(2n-1)}~{},$
(29)

where $z=\rho\sigma(v-v_{0})^{2}\exp{(\delta+\beta v)},$ with
$\gamma=a+\beta,\ \delta=\alpha-\beta v_{0},\ \rho=(a^{2}+6a\beta+\beta^{2}),$
and $\sigma=\gamma/\beta$.

A ratio test of the series (29) shows that it is convergent for any value of
$z$. The summation can be approximated as

$\sum_{n=0}^{\infty}\frac{z^{n}}{n!(2n-1)}=\frac{4}{z^{2}}(e^{z}-1)-\frac{4}{z}-3+\frac{z}{3}~{},$
(30)

which is an excellent approximation because the series expansion terms on the
right hand side $(b_{n}z^{n})$ differ from those on the left hand side
$(a_{n}z^{n})$ within $|b_{n}-a_{n}|\sim 10^{-n-3}$.

Thus the mass function of the black hole, as measured by an infalling
observer, is obtained from equations (29) and equation (30) as

$\begin{split}&\frac{m(v)-m_{0}}{r_{0}^{2}}\\\ &\hskip
22.76219pt\approx\frac{4\exp{\left[(v-v_{0})^{2}\rho\sigma e^{(\delta+\gamma
v)}-2(\delta+\gamma v)\right]}}{(v-v_{0})^{5}\rho^{2}\sigma^{2}}\\\ &\hskip
42.67912pt+\frac{1}{3}(v-v_{0})\rho e^{(\delta+\gamma v)}\\\ \end{split}$ (31)

where only the dominant contributions have been included.

The first term on the right hand side of equation (31) gives the most dominant
contribution because of its double exponential nature. Consequently, the CH
instability is enormously enhanced in the limit $v\xrightarrow[]{}\infty$.

Figure 2: Plot displaying mass superinflation with different values of scalar
mass $\mu$.

##  6 Discussion and Conclusion

In this work, we considered the dynamics of a massive scalar field coupled to
the dynamics of the mass function of the Reissner-Nordström black hole to find
the nature of Cauchy horizon (CH) singularity. The scalar field $\Phi(r,v)$
was assumed to be chargeless and spherically symmetric with its intensity
falling off like $1/r^{2}$. As a result of this perturbation, the redshift
function modifies to $f=1-2[m_{0}+m(v)]/r+Q^{2}/r^{2}$. Employing the
Einstein-Maxwell field equation coupled with the Klein-Gordon equation, we
obtained the coupled dynamics of the mass function $m(v)$ and the scalar field
$\Phi(r,v)$ given by equations (6) and (7). To solve the coupled set of
nonlinear differential equations, we developed a perturbative method and
obtained a second order nonlinear differential equation given by equation
(20). To solve this equation, we employed the Adomian decomposition method
(ADM). This results in a double exponential behaviour for the mass function
$m(v)$ in the asymptotic limit, as shown in equation (31).

It is evident from equation (31) that the mass of the Reissner-Nordström black
hole, as measured by an infalling observer, shows a rapid and unbounded double
exponential growth upon approaching the CH.

Figure 2 displays the behaviour of the mass function for different values of
scalar mass $\mu$ including the massless case. The plot illustrates that the
mass inflation in the massive case is stronger than in the massless case, the
behaviour being dictated by the double exponential nature in both cases. This
is also evident from the expression for the mass function $m(v)$ given by
equation (31) [where the coefficients are $\rho=(a^{2}+6a\beta+\beta^{2})$ and
$\gamma=a+\beta$, with $a=\mu^{2}r_{0}^{2}$, $\sigma=\gamma/\beta$,
$\beta=m_{0}/r_{0}^{2}$ and $\delta=\alpha-\beta v_{0}$ is a constant]. It is
important to note that such double exponential growth is a enormously stronger
mass inflation compared to those obtained previous studies.

In the near horizon and long time limits ($x\xrightarrow[]{}0$ and
$v\xrightarrow[]{}\infty$), the dominant behaviour of the mass function and
the scalar field can be expressed as

$m(v)\sim v^{-5}e^{-2\gamma v}\exp{[v^{2}\rho\sigma e^{\gamma v}]}$

and

$\phi(v)\sim\exp\bigg{[}\frac{1}{2r_{0}^{2}}\int m(v)\,dv\bigg{]}.$

The strongly diverging nature of the scalar field at the CH comes from the
strong gravitational blueshift of the continuous scalar field surrounding the
black hole. Because of this blueshift, the scalar field rapidly accumulates at
the CH with increasing energy which, in turn, results in a rapid double
exponential growth in the mass of the black hole, creating what may be termed
as mass superinflation. This implies that the information accumulated at the
inner horizon would be seen by the infalling observer and all the history of
the external universe containing the scalar field will be available to him/her
within a fraction of time, as was first predicted by Penrose [6]. Because of
the double exponentially diverging mass function, the Kretschmann scalar
$R_{\mu\nu\alpha\beta}R^{\mu\nu\alpha\beta}\propto m^{2}$ is also a double
exponentially diverging function, confirming an enormously strong singularity
at the CH. Such strong divergence suggests that there is no possibility of
extending the spacetime beyond the Cauchy horizon.

##  7 Acknowledgement

Nihar Ranjan Ghosh is supported by a Research Fellowship of the Ministry of
Education (MoE), Government of India.

## References

  * [1] Hans Reissner.  Über die eigengravitation des elektrischen feldes nach der einsteinschen theorie.  Annalen der Physik, 355(9):106–120, 1916.
  * [2] G Nordstrom.  On the energy of the gravitational field in einstein’s theory ii verhandl.  Koninkl. Ned. Akad. Wetenschap., Afdel. Natuurk., Amsterdam, 26:1201–1208, 1918.
  * [3] Roger Penrose.  Gravitational collapse: The role of general relativity.  Nuovo Cimento Rivista Serie, 1:252, 1969.
  * [4] ID Novikov.  Change of relativistic collapse into anticollapse and kinematics of a charged sphere.  JETP Lett.(USSR)(Engl. Transl.), 3, 1966.
  * [5] Michael Simpson and Roger Penrose.  Internal instability in a reissner-nordström black hole.  International Journal of Theoretical Physics, 7:183–197, 1973.
  * [6] R. Penrose.  Battelle rencontres ‐ 1967 lectures in mathematics and physics (ed. by c. m. dewitt and j. a. wheeler). xvii + 557 s. m. fig. new york/amsterdam 1968\. w. a. benjamin, inc.  Zamm-zeitschrift Fur Angewandte Mathematik Und Mechanik, 50:782–782, 1970.
  * [7] Amos Ori.  Inner structure of a charged black hole: An exact mass-inflation solution.  Physical review letters, 67(7):789, 1991.
  * [8] Andrew JS Hamilton and Pedro P Avelino.  The physics of the relativistic counter-streaming instability that drives mass inflation inside black holes.  Physics reports, 495(1):1–32, 2010.
  * [9] Eric Poisson and Werner Israel.  Inner-horizon instability and mass inflation in black holes.  Physical review letters, 63(16):1663, 1989.
  * [10] Werner Israel and Eric Poisson.  Internal structure of black holes.  Physical Review D, 41(6):1796, 1990.
  * [11] Eric G Brown, Robert Mann, and Leonardo Modesto.  Mass inflation in the loop black hole.  Physical Review D, 84(10):104041, 2011.
  * [12] Richard H Price.  Nonspherical perturbations of relativistic gravitational collapse. i. scalar and gravitational perturbations.  Physical Review D, 5(10):2419, 1972.
  * [13] Frank J Tipler.  Singularities in conformally flat spacetimes.  Physics Letters A, 64(1):8–10, 1977.
  * [14] George FR Ellis and Bernd G Schmidt.  Singular space-times.  General Relativity and Gravitation, 8:915–953, 1977.
  * [15] Lior M Burko and Amos Ori.  Analytic study of the null singularity inside spherical charged black holes.  Physical Review D, 57(12):R7084, 1998.
  * [16] M Dafermos and J Luk.  The interior of dynamical vacuum black holes i: the, 2017.
  * [17] Paul CW Davies, Stephen A Fulling, and William G Unruh.  Energy-momentum tensor near an evaporating black hole.  Physical Review D, 13(10):2720, 1976.
  * [18] William A Hiscock.  Stress-energy tensor near a charged, rotating, evaporating black hole.  Physical Review D, 15(10):3054, 1977.
  * [19] Roberto Balbinot and Eric Poisson.  Mass inflation: The semiclassical regime.  Physical review letters, 70(1):13, 1993.
  * [20] Nicholas David Birrell and Paul Charles William Davies.  Quantum fields in curved space.  Cambridge university press, 1984.
  * [21] Dong-il Hwang and Dong-han Yeom.  Internal structure of charged black holes.  Physical Review D, 84(6):064020, 2011.
  * [22] Carlos Barceló, Valentin Boyanov, Raúl Carballo-Rubio, and Luis J Garay.  Black hole inner horizon evaporation in semiclassical gravity.  Classical and Quantum Gravity, 38(12):125003, 2021.
  * [23] Carlos Barceló, Valentin Boyanov, Raúl Carballo-Rubio, and L J Garay.  Classical mass inflation versus semiclassical inner horizon inflation.  Physical Review D, 106(12):124006, 2022.
  * [24] Raúl Carballo-Rubio, Francesco Di Filippo, Stefano Liberati, Costantino Pacilio, and Matt Visser.  Regular black holes without mass inflation instability.  Journal of High Energy Physics, 2022(9):1–14, 2022.
  * [25] Edgardo Franzin, Stefano Liberati, Jacopo Mazza, and Vania Vellucci.  Stable rotating regular black holes.  Physical Review D, 106(10):104060, 2022.
  * [26] Tyler McMaken.  Semiclassical instability of inner-extremal regular black holes.  Physical Review D, 107(12):125023, 2023.
  * [27] PR Brady and E Poisson.  Cauchy horizon instability for reissner-nordstrom black holes in de sitter space.  Classical and Quantum Gravity, 9(1):121, 1992.
  * [28] Felicity Mellor and Ian Moss.  Stability of black holes in de sitter space.  Physical Review D, 41(2):403, 1990.
  * [29] Patrick R Brady, Dario Nunez, and Sukanya Sinha.  Cauchy horizon singularity without mass inflation.  Physical Review D, 47(10):4239, 1993.
  * [30] Rong-Gen Cai and Ru-Keng Su.  Black holes in de sitter space and the stability conjecture of cauchy horizons.  Physical Review D, 52(2):666, 1995.
  * [31] TM Helliwell and DA Konkowski.  Testing a stability conjecture for cauchy horizons.  Physical Review D, 47(10):4322, 1993.
  * [32] DA Konkowski and TM Helliwell.  Instabilities of the cauchy horizon in kerr black holes.  Physical Review D, 50(2):841, 1994.
  * [33] Dragoljub Marković and Eric Poisson.  Classical stability and quantum instability of black-hole cauchy horizons.  Physical review letters, 74(8):1280, 1995.
  * [34] G Adomian.  A review of the decomposition method and some recent results for nonlinear equations.  Mathematical and Computer Modelling, 13(7):17–43, 1990.
  * [35] George Adomian.  A review of the decomposition method in applied mathematics.  Journal of mathematical analysis and applications, 135(2):501–544, 1988.
  * [36] George Adomian.  Solving frontier problems of physics: the decomposition method, volume 60.  Springer Science & Business Media, 2013.
  * [37] George Adomian.  Nonlinear stochastic operator equations.  Academic press, 2014.
  * [38] Jun-Sheng Duan, Randolph Rach, Dumitru Baleanu, and Abdul-Majid Wazwaz.  A review of the adomian decomposition method and its applications to fractional differential equations.  Communications in Fractional Calculus, 3(2):73–99, 2012.