# Fully Hyperbolic Graph Convolution Network for Recommendation

Liping Wang1,2,*, Fenyu Hu1,2*, Shu Wu1,2,†, and Liang Wang1,2 1Center for
Research on Intelligent Perception and Computing, Institute of Automation,
Chinese Academy of Sciences2School of Artificial Intelligence, University of
Chinese Academy of Sciences<EMAIL_ADDRESS><EMAIL_ADDRESS>shu.wu<EMAIL_ADDRESS>
(2018)

###### Abstract.

Recently, Graph Convolution Network (GCN) based methods have achieved
outstanding performance for recommendation. These methods embed users and
items in Euclidean space, and perform graph convolution on user-item
interaction graphs. However, real-world datasets usually exhibit tree-like
hierarchical structures, which make Euclidean space less effective in
capturing user-item relationship. In contrast, hyperbolic space, as a
continuous analogue of a tree-graph, provides a promising alternative. In this
paper, we propose a fully hyperbolic GCN model for recommendation, where all
operations are performed in hyperbolic space. Utilizing the advantage of
hyperbolic space, our method is able to embed users/items with less distortion
and capture user-item interaction relationship more accurately. Extensive
experiments on public benchmark datasets show that our method outperforms both
Euclidean and hyperbolic counterparts and requires far lower embedding
dimensionality to achieve comparable performance.

recommendation system, graph neural networks, hyperbolic space

††copyright: acmcopyright††journalyear: 2018††doi:
10.1145/1122445.1122456††conference: Woodstock ’18: ACM Symposium on Neural
Gaze Detection; June 03–05, 2018; Woodstock, NY††booktitle: Woodstock ’18: ACM
Symposium on Neural Gaze Detection, June 03–05, 2018, Woodstock, NY††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs: Computer systems organization
Embedded systems††ccs: Computer systems organization Redundancy††ccs: Computer
systems organization Robotics††ccs: Networks Network reliability

##  1\. Introduction

In information era, recommendation systems have been widely adopted to perform
personalized information filtering (Covington et al., 2016; Ying et al.,
2018). Even though there are many recommendation paradigms, collaborative
filtering (Ebesu et al., 2018; Wang et al., 2019) which generates
recommendations by utilizing available historical interactions, remains a
fundamental and challenging task. The core idea behind collaborative filtering
is to learn compact user/item embeddings and infer a user’s preference to an
item according to the distance between their embeddings.

From the perspective of graph, user-item interactions can be viewed as a bi-
partite graph (Bollobás, 2013), where nodes represent users/items and edges
represent their interactions. As a powerful tool of analyzing graph-structured
data, Graph Neural Networks(GNNs) (Kipf and Welling, 2017; Veličković et al.,
2018; Hamilton et al., 2017) have recently demonstrated great success across
various domains, including recommendation systems. Employing multiple layers
of neighborhood aggregation, GNN-based methods (Wang et al., 2019; He et al.,
2020) have achieved the state-of-the-art performance on diverse public
benchmarks.

Although GNN-based methods have achieved outstanding performance, it might not
be appropriate to adopt Euclidean space to embed users and items. In real-
world scenarios, user-item bi-partite graphs usually exhibit tree-like
hierarchical structures (Adcock et al., 2013), in which the number of a node’s
neighbors grows exponentially with respect to the number of hops. Ideally,
neighbors of node $v$ should be embedded in the ball centering at $v$’s
embedding, and the distance between embeddings should reflect the number of
hops between nodes. Nevertheless, in Euclidean space, the volume of a ball
only grows polynomially as a function of radius. Hence, embedding
exponentially-growing number of neighbors into polynomially-growing size of
volume would make distance between embeddings less accurate to reflect
distance between nodes in the graph, and this is called distortion (Chami et
al., 2019). This kind of distortion makes it difficult to infer a user’s
preference to a target item according to the distance between their
embeddings.

In contrast, hyperbolic space (Chami et al., 2019) in which the volume of a
ball grows exponentially with radius offers a promising alternative. Compared
with Euclidean space, hyperbolic space is more suitable for modeling user-item
interaction graphs which exhibit strong tree-like hierarchical structures.
Accordingly, it is a natural choice to conduct user/item embedding and graph
convolution in hyperbolic space for recommendation. To the best of our
knowledge, the only work adopting a similar idea is HGCF (Sun et al., 2021).
However, resorting to tangent space to realize graph convolution makes the
performance of HGCF inferior for the two following reasons. On the one hand,
tangent space is only local linear approximation of hyperbolic space. During
the process of message propagation, errors caused by approximation accumulate
and spread to the whole graph. As a result, influence from high-order
neighbors cannot be captured accurately. On the other hand, tangent space is
actually Euclidean space according to its definition (Boothby, 1986). Hence,
the advantage of hyperbolic space in modeling user-item interaction
relationship with less distortion cannot be fully utilized.

To overcome the limitation of Euclidean space and obtain more accurate
user/item embeddings, we design a novel fully hyperbolic GCN framework
specially for collaborative filtering. All operations are conducted in
hyperbolic space, more specifically, in the Lorentz model, and we name it
Lorentz Graph Collaborative Filtering (LGCF).

The main contributions of this work are summarized as follows:

  * •

We propose a fully hyperbolic graph convolution network for recommendation.

  * •

We conduct extensive experiments on multiple public benchmark datasets, and
the results demonstrate the superiority of our method.

##  2\. Preliminaries

Problem Formulation. In this paper, the standard collaborative filtering set-
up is considered. Let $U=\\{u_{1},u_{2},\dots,u_{n}\\}$ be the set of $n$
users, and $I=\\{i_{1},i_{2},\dots,i_{m}\\}$ be the set of $m$ items.
Historical interactions between users and items are represented as a binary
matrix $B$, where $B_{i,j}=1$ if the $i$-th user interacts with the $j$-th
item, otherwise $0$. Given historical interactions $B$, the goal is to predict
potential interactions.

GNN-based Recommendation Methods. GNN-based recommendation methods have
received increasing attention for their ability to learn rich node
representations. In collaborative filtering setting, GCN has replaced matrix
factorization (Rendle et al., 2009) and shown leading performance. In order to
capture the influence from high-order neighbors, Wang et al. (Wang et al.,
2019) proposed NGCF, a GCN framework performing on user-item interaction bi-
partite graphs. Subsequently, He et al. (He et al., 2020) empirically found
that two most common designs in GCNs — feature transformation and nonlinear
activation, contribute little to the performance of collaborative filtering.
Hence, they proposed a simplified architecture – LightGCN (He et al., 2020).

##  3\. Our Method

Figure 1. Illustration of LGCF. First, users and items are embedded into the
Lorentz model of hyperbolic space. Then, multiple graph convolution layers
(only one shown in the figure for simplicity) are adopted to aggregate
information from neighbors. In each layer, embeddings are first mapped by $F$
to the Klein model in which graph convolution is performed. After that,
$F^{-1}$ maps refined embeddings back to the Lorentz model. Finally, LGCF
infers a user’s preference to an item according to the distance between their
embeddings.

As illustrated in Figure 1, there are three components in LGCF: (1) an
embedding layer that provides and initializes user/item embeddings in
hyperbolic space; (2) multiple graph convolution layers that propagate
user/item embeddings over the graph; and (3) a prediction layer that estimates
a user’s preference to an item by computing the distance between their
embeddings.

###  3.1. Embed Users/Items in Hyperbolic Space

Existing GNN-based recommendation methods usually embed users and items in the
same Euclidean space. To accurately model user-item interaction relationship,
we investigate the utilization of hyperbolic space. There are several models
of hyperbolic space, such as the Lorentz model, the Klein model and the
Poincaré ball model. In this paper, we choose the Lorentz model due to its
simplicity and numerical stability. Formally, the Lorentz model of
$d$-dimensional hyperbolic space is defined as:

(1)
$\mathcal{L}=\\{\boldsymbol{x}=[x_{0},x_{1},\cdots,x_{d}]\in\mathbb{R}^{d+1}:\left<\boldsymbol{x},\boldsymbol{x}\right>_{\mathcal{L}}=-1,x_{0}>0\\},$

where $\left<\boldsymbol{x},\boldsymbol{y}\right>_{\mathcal{L}}$ is Lorentz
inner product and is defined as
$\left<\boldsymbol{x},\boldsymbol{y}\right>_{\mathcal{L}}=-x_{0}y_{0}+\sum_{i=1}^{d}x_{i}y_{i}$.
In addition, at an arbitrary point $\boldsymbol{x}\in\mathcal{L}$, hyperbolic
space can be locally approximated by a linear Euclidean space. And this
approximated Euclidean space is termed as tangent space
$\mathcal{T}_{\boldsymbol{x}}\mathcal{L}$ in which norm
$\left\lVert\boldsymbol{x}\right\rVert_{\mathcal{L}}=\left<\boldsymbol{x},\boldsymbol{y}\right>$
is well defined. In LGCF, both users and items are embedded in the same
Lorentz model of hyperbolic space.

It is well known that random initialization can have a significant impact on
optimization in training (Sutskever et al., 2013). A common practice in
Euclidean space is Gaussian distribution initialization. Similarly, we design
an initialization strategy for embeddings based on Wrapped Normal Distribution
(Nagano et al., 2019) which is a generalization of Gaussian distribution to
hyperbolic space.

###  3.2. Graph Convolution Layer

The basic idea of GCN-based recommendation models is learning representations
for users and items by aggregating neighbors’ information iteratively over the
interaction graph. In order to apply GCN for recommendation where users and
items are embedded in hyperbolic space, we design graph convolution layers
specially since naïve generalization will drive embeddings out of hyperbolic
space. Before that, we give a brief review of existing graph convolution
layers in Euclidean space.

In Euclidean space, a graph convolution operation is composed of three steps:
feature transformation, neighborhood aggregation and non-linear activation.
Among them, feature transformation is performed through linear transformation.
In LGCF, we discard the feature transformation operation for two reasons. On
the one hand, different from attributed graphs (e.g., citation networks) where
nodes bring rich feature information, nodes in user-item interaction graphs
contain no semantics but one-hot IDs. In this case, feature transformation may
provide no benefits, and could bring difficulties to training. On the other
hand, linear transformation (matrix-vector multiplication) is not well-defined
in hyperbolic space since it is not a vector space.

Neighborhood Aggregation. Exising neighborhood aggregation of the $l$-th layer
can be summarized as:

(2)
$\boldsymbol{e}^{(l)}_{i}=\sum_{j\in\hat{\mathcal{N}}(i)}w_{ij}\boldsymbol{e}_{j}^{(l-1)},$

in which $\boldsymbol{e}_{j}^{(l-1)}$ represents the $j$-th node’s embedding
and $\hat{\mathcal{N}}(i)$ denotes the set consisting of node $i$ and its
neighborhood nodes. A natural generalization of mean to hyperbolic space is
Einstein midpoint (Ungar, 2008). However, it is defined in the Klein model
$\mathcal{K}=\left\\{\boldsymbol{x}\in\mathbb{R}^{d}:\|\boldsymbol{x}\|<1\right\\}$,
while user/item embeddings lie in the Lorentz model. Luckily, there are
isometric maps between the Klein model and the Lorentz model defined as
follows:

(3) $\displaystyle\begin{split}F(\boldsymbol{x})&=\frac{[x_{1},x_{2},\dots
x_{d}]}{x_{0}},\\\
F^{-1}(\boldsymbol{k})&=\frac{1}{1-\left\lVert\boldsymbol{k}\right\rVert^{2}}[1,\boldsymbol{k}],\end{split}$

in which $\boldsymbol{x}=[x_{0},x_{1},\dots,x_{d}]\in\mathcal{L}$ and
$\boldsymbol{k}\in\mathcal{K}$.

Hence, we propose a neighborhood aggregation strategy utilizing the Klein
model as an intermediate bridge. Specifically, neighborhood aggregation can be
divided into three steps. First, current embeddings
$\\{\boldsymbol{e}_{1},\boldsymbol{e}_{2},\dots,\boldsymbol{e}_{n}\\}$ are
mapped from the Lorentz model to the Klein model by
$\boldsymbol{k}_{i}=F(\boldsymbol{e}_{i})$. Then, neighborhood aggregation is
conducted in the Klein model as follows:

(4)
$\displaystyle\begin{split}\gamma_{i}=\frac{1}{\sqrt{1-\left\lVert\boldsymbol{k}_{j}\right\rVert^{2}}},\\\
\boldsymbol{k}_{i}^{\prime}=\frac{\sum_{j\in\hat{N}(i)}\gamma_{j}\boldsymbol{k}_{j}}{\sum_{j\in\hat{N}(i)}\gamma_{j}}.\\\
\end{split}$

Here, we obtain aggregated user/item embeddings in the Klein model. Last,
aggregated embeddings $\boldsymbol{k}_{i}^{\prime}$ are mapped back to the
Lorentz model through
$\boldsymbol{z}_{i}=F^{-1}(\boldsymbol{k}_{i}^{\prime})$.

Nonlinear Activation Layer. In Euclidean space, nonlinear activation has
proven to be a key component in modern neural networks. However, direct
adoption will drive the computation result out of hyperbolic space. To fix
this problem, we design a calibration strategy following general activation.
Formally, let
$\boldsymbol{x_{i}}=[x_{0},x_{1},\dots,x_{d}]=\sigma(\boldsymbol{z}_{i})$ be
the output of general activation function $\sigma$, e.g., ReLU. The first
element of $\boldsymbol{x}$ is calibrated while the other elements remain
unchanged to pull the activated embedding back to hyperbolic space:

(5)
$\boldsymbol{e}_{i}^{\prime}=f_{c}(\boldsymbol{x}_{i})=\left[\sqrt{1+\sum_{j=1}^{d}x_{j}^{2}},x_{1},x_{2},\dots,x_{d}\right].$

###  3.3. Prediction Layer

After propagating with $L$ graph convolution layers, we obtain multiple
representations for users and items. Representations generated in different
layers emphasize messages passed through different connections, and they
reflect users’ preference or items’ attributes from different perspectives.
Recommendation models operating in hyperbolic space usually estimate a user’s
preference to a target item according to the distance or similarity metric
between their representations. In the Lorentz model, the generalization of a
straight line in Euclidean space is geodesics which gives the shortest
distance between two points $\boldsymbol{x},\boldsymbol{y}$. Formally, the
geodesics distance between $\boldsymbol{x}$ and $\boldsymbol{y}$ is defined
as:

(6)
$d_{\mathcal{L}}(\boldsymbol{x},\boldsymbol{y})=\operatorname{arcosh}(-\left<\boldsymbol{x},\boldsymbol{y}\right>_{\mathcal{L}}).$

Hence, in LGCF, we infer the preference of user $u$ to item $i$ based on the
geodesics distance between their corresponding representations. Further, to
utilize different semantics captured by different layers, we take
representations learned by different layers into consideration simultaneously.
In summary, LGCF estimates the preference of user $u$ to target item $i$ as:

(7)
$y(u,i)=\frac{1}{\sum_{l=1}^{L}d_{\mathcal{L}}^{2}(\boldsymbol{e}_{u}^{(l)},\boldsymbol{e}_{i}^{(l)})},$

in which $\boldsymbol{e}_{u}^{(l)}$ and $\boldsymbol{e}_{i}^{(l)}$ are
representations generated by the $l$-th layer. Even though there are multiple
choices for layer aggregation, such as weighted average, max pooling, LSTM,
etc., we find simple summation adopted here works well empirically.

###  3.4. Margin Ranking Loss

Margin ranking loss (Tay et al., 2018) has been a competitive choice for
distance-based recommendation systems. Since it encourages positive and
negative user-item pairs to be separated up to a given margin. Once the
difference between a negative user-item pair and a positive one is greater
than the margin, these two user-item pairs will make no contribution to the
loss. In this way, hard pairs violating the margin are focused all the time,
making optimization much easier. We extend margin ranking loss to hyperbolic
space based on geodesics distance. Given a sampled positive user-item pair
$(u,i)$ and a negative one $(u,j)$, geodesics margin loss is defined as:

(8)
$\ell_{g}(u,i,j)=\sum_{l=1}^{L}\min\left(d_{\mathcal{L}}^{2}(\boldsymbol{e}_{u}^{(l)},\boldsymbol{e}_{j}^{(l)})-d_{\mathcal{L}}^{2}(\boldsymbol{e}_{u}^{(l)},\boldsymbol{e}_{i}^{(l)})-m,0\right),$

where $m$ is a non-negative hyper-parameter. Note that representations
obtained by different layers contribute to the loss simultaneously. This not
only makes it possible to utilize different semantics captured by different
layers, but also decreases the difficulty of optimization due to the residual
connection (He et al., 2016).

###  3.5. Optimization

The only parameter of LGCF is the embedding matrix of users and items. These
embeddings lie in the Lorentz model of hyperbolic space which is out of the
range of common optimization algorithms such as SGD (Robbins and Monro, 1951)
and Adam (Kingma and Ba, 2014). Hence, we employ RGSD (Bonnabel, 2013), a
generation of SGD to hyperbolic space, which mimics SGD’s behavior while
taking into account the geometry of hyperbolic space.

##  4\. Experiments

###  4.1. Set Up

Table 1. Dataset statistics. Dataset | #User | #Item | #Interactions
---|---|---|---
Amazon-CD | 22,947 | 18,395 | 422,301
Amazon-Book | 52,406 | 41,264 | 1,861,118
Yelp2020 | 91,174 | 45,063 | 1,940,014

Datasets and Baselines. Following HGCF (Sun et al., 2021), we employ Amazon-CD
(Ni et al., 2019), Amazon-Book (Ni et al., 2019) and Yelp2020 (yel, 2020)
datasets. Dataset statistics are provided in Table 1. Each dataset is split
into 80-20 train and test sets. Multiple competitive baseline methods from
three categories are compared: BPRMF (Rendle et al., 2009), NGCF (Wang et al.,
2019), LightGCN (He et al., 2020), HVAE (Mirvakhabova et al., 2020), HGCF (Sun
et al., 2021). Among them, BPRMF optimizes matrix factorization by Bayesian
personalize ranking (BPR) loss (Rendle et al., 2009). NGCF and LightGCN employ
GNN in Euclidean space. HVAE combines variational auto-encoder(VAE) with
hyperbolic geometry. Last, HGCF applies the latest Hyperbolic GCN (Chami et
al., 2019) to recommendation systems.

Implementation. For a fair comparison, the embedding dimensionality is set to
50 for all methods, and the same negative sampling strategy is adopted. For
all baseline methods, suggested settings in original papers are followed. Grid
search for hyper-parameters are conducted following HGCF (Sun et al., 2021).
For LGCF, the number of GCN layers is set to 3. We set the learning rate to
0.001 and weight decay to 0.005. And model is trained for 1000 epochs. Hyper-
parameter margin $m$ is tuned from $[0.01,2]$.

###  4.2. Overall Results

Table 2. Recall results for all datasets.

Datasets | Metrics | BPRMF | NGCF | LightGCN | HVAE | HGCF | LGCF
---|---|---|---|---|---|---|---
Amazon-CD | R@10 | 0.0779 | 0.0758 | 0.0929 | 0.0781 | 0.0962 | 0.0996
R@20 | 0.1200 | 0.1150 | 0.1404 | 0.1147 | 0.1455 | 0.1503
Amazon-Book | R@10 | 0.0611 | 0.0658 | 0.0799 | 0.0774 | 0.0867 | 0.0899
R@20 | 0.0794 | 0.1050 | 0.1248 | 0.1125 | 0.1318 | 0.1360
Yelp2020 | R@10 | 0.0325 | 0.0458 | 0.0522 | 0.0421 | 0.0543 | 0.0573
R@20 | 0.0556 | 0.0764 | 0.0866 | 0.0691 | 0.0884 | 0.0946

Table 3. NDCG results for all datasets.

Datasets | Metrics | BPRMF | NGCF | LightGCN | HVAE | HGCF | LGCF
---|---|---|---|---|---|---|---
Amazon-CD | N@10 | 0.0610 | 0.0591 | 0.0726 | 0.0629 | 0.0751 | 0.0780
N@20 | 0.0974 | 0.0718 | 0.0881 | 0.0749 | 0.0909 | 0.0945
Amazon-Book | N@10 | 0.0594 | 0.0655 | 0.0780 | 0.0778 | 0.0869 | 0.0906
N@20 | 0.0971 | 0.0791 | 0.0938 | 0.0901 | 0.1022 | 0.1063
Yelp2020 | N@10 | 0.0283 | 0.0405 | 0.0461 | 0.0371 | 0.0458 | 0.0485
N@20 | 0.0512 | 0.0513 | 0.0582 | 0.0465 | 0.0585 | 0.0612

Recall and NDCG results for all datasets are reported in Table 2 and Table 3
respectively. We can see that LGCF consistently outperforms other methods on
all the three datasets. Compared with LightGCN, LGCF achieves an improvement
up to 16.15% in NDCG@10 on Amazon-Book dataset, demonstrating the superiority
of hyperbolic space over Euclidean space in modeling real-world user-item
interactions.

Among the baseline methods, HGCF is the most competitive counterpart of LGCF.
Even though HGCF adopts hyperbolic space, it resorts to tangent space to
conduct aggregation operations, which brings inevitable distortion. In
contrast, LGCF performs all graph convolution operations in hyperbolic space.
LGCF outperforms HGCF with wide margins on all the datasets, showing the
information loss introduced by tangent space in HGCF.

###  4.3. Ablation Study

Table 4. Results for LGCF and LGCF-tangent on Yelp2020. | R@10 | R@20 | N@10 | N@20
---|---|---|---|---
LGCF | 0.0573 | 0.0946 | 0.0485 | 0.0612
LGCF-tangent | 0.0545 | 0.0895 | 0.0463 | 0.0586

To further analyze the effect of fully hyperbolic graph convolution network,
we conduct an ablation study on Yelp2020, the largest one among the three
datasets. Since simply replacing fully hyperbolic graph convolution with
regular Euclidean graph convolution would drive user/item embeddings out of
hyperbolic space, we conduct convolution operations in tangent space, and name
this model variant as LGCF-tangent. From experimental results shown in Table
4, we can observe that there is a wide margin between the performance of LGCF
and LGCF-tangent. This is because, in LGCF-tangent, errors caused by tangent
space approximation accumulate and spread to the whole graph. As a result,
influence from neighbors, especially high-order neighbors, cannot be captured
accurately.

###  4.4. Effect of Embedding Dimensionality

Figure 2. Recall@20 and NDCG@20 on Amazon-CD dataset with an embedding
dimensionality varying from 20 to 50.

In order to validate the advantage of hyperbolic space to learn compact
representations, we compare the results of LGCF and LightGCN with different
values of embedding dimensionality. From Figure 2, we can observe that LGCF
outperforms LightGCN consistently at all dimensionality values, and the
greatest margin occurs at lower dimensionality. LGCF requires far lower
embedding dimensionality to achieve comparable performance to its Euclidean
analogue. This reflects that LGCF’s advantage is more prominent when the
embedding dimensionality cannot be large due to limited computing and storage
resources.

##  5\. Conclusion

In this paper, we propose LGCF, a fully hyperbolic GCN model for
recommendation. Utilizing the advantage of hyperbolic space, LGCF is able to
embed users/items with less distortion and capture user-item interaction
relationship more accurately. Extensive experiments on public benchmark
datasets show that LGCF outperforms both Euclidean and hyperbolic counterparts
and requires far lower embedding dimensionality to achieve comparable
performance.

## References

  * (1)
  * yel (2020) 2020\.  _Yelp Open Dataset_.  https://www.yelp.com/dataset
  * Adcock et al. (2013) Aaron B Adcock, Blair D Sullivan, and Michael W Mahoney. 2013\.  Tree-like structure in large social and information networks. In _IEEE International Conference on Data Mining_.
  * Bollobás (2013) Béla Bollobás. 2013\.  _Modern graph theory_.
  * Bonnabel (2013) Silvere Bonnabel. 2013\.  Stochastic gradient descent on Riemannian manifolds.  _IEEE Trans. Automat. Control_ (2013).
  * Boothby (1986) William M Boothby. 1986\.  _An introduction to differentiable manifolds and Riemannian geometry_.
  * Chami et al. (2019) Ines Chami, Rex Ying, Christopher Ré, and Jure Leskovec. 2019\.  Hyperbolic graph convolutional neural networks.  _Advances in Neural Information Processing Systems_ (2019).
  * Covington et al. (2016) Paul Covington, Jay Adams, and Emre Sargin. 2016\.  Deep neural networks for youtube recommendations. In _Proceedings of the 10th ACM Conference on Recommender Systems_.
  * Ebesu et al. (2018) Travis Ebesu, Bin Shen, and Yi Fang. 2018.  Collaborative memory network for recommendation systems. In _The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval_.
  * Hamilton et al. (2017) William L. Hamilton, Rex Ying, and Jure Leskovec. 2017\.  Inductive Representation Learning on Large Graphs. In _Advances in Neural Information Processing Systems_.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. 2016\.  Deep Residual Learning for Image Recognition. In _IEEE Conference on Computer Vision and Pattern Recognition_.
  * He et al. (2020) Xiangnan He, Kuan Deng, Xiang Wang, Yan Li, Yongdong Zhang, and Meng Wang. 2020\.  Lightgcn: Simplifying and powering graph convolution network for recommendation. In _Proceedings of the 43rd International ACM SIGIR Conference on Research and Development in Information Retrieval_.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * Kipf and Welling (2017) Thomas N. Kipf and Max Welling. 2017.  Semi-Supervised Classification with Graph Convolutional Networks. In _International Conference on Learning Representations_.
  * Mirvakhabova et al. (2020) Leyla Mirvakhabova, Evgeny Frolov, Valentin Khrulkov, Ivan Oseledets, and Alexander Tuzhilin. 2020\.  Performance of Hyperbolic Geometry Models on Top-N Recommendation Tasks. In _Fourteenth ACM Conference on Recommender Systems_.
  * Nagano et al. (2019) Yoshihiro Nagano, Shoichiro Yamaguchi, Yasuhiro Fujita, and Masanori Koyama. 2019.  A wrapped normal distribution on hyperbolic space for gradient-based learning. In _International Conference on Machine Learning_.
  * Ni et al. (2019) Jianmo Ni, Jiacheng Li, and Julian McAuley. 2019.  Justifying recommendations using distantly-labeled reviews and fine-grained aspects. In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing_.
  * Rendle et al. (2009) Steffen Rendle, Christoph Freudenthaler, Zeno Gantner, and Lars Schmidt-Thieme. 2009.  BPR: Bayesian Personalized Ranking from Implicit Feedback. In _Proceedings of the Twenty-Fifth Conference on Uncertainty in Artificial Intelligence_.
  * Robbins and Monro (1951) Herbert Robbins and Sutton Monro. 1951.  A stochastic approximation method.  _The Annals of Mathematical Statistics_ (1951).
  * Sun et al. (2021) Jianing Sun, Zhaoyue Cheng, Saba Zuberi, Felipe Pérez, and Maksims Volkovs. 2021.  HGCF: Hyperbolic Graph Convolution Networks for Collaborative Filtering.  (2021).
  * Sutskever et al. (2013) Ilya Sutskever, James Martens, George Dahl, and Geoffrey Hinton. 2013.  On the importance of initialization and momentum in deep learning. In _International Conference on Machine Learning_.
  * Tay et al. (2018) Yi Tay, Luu Anh Tuan, and Siu Cheung Hui. 2018.  Latent relational metric learning via memory-based attention for collaborative ranking. In _Proceedings of the 2018 World Wide Web Conference_.
  * Ungar (2008) Abraham Albert Ungar. 2008\.  A gyrovector space approach to hyperbolic geometry.  _Synthesis Lectures on Mathematics and Statistics_ (2008).
  * Veličković et al. (2018) Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. 2018.  Graph Attention Networks.  _International Conference on Learning Representations_ (2018).
  * Wang et al. (2019) Xiang Wang, Xiangnan He, Meng Wang, Fuli Feng, and Tat-Seng Chua. 2019.  Neural graph collaborative filtering. In _Proceedings of the 42nd international ACM SIGIR Conference on Research and Development in Information Retrieval_.
  * Ying et al. (2018) Rex Ying, Ruining He, Kaifeng Chen, Pong Eksombatchai, William L Hamilton, and Jure Leskovec. 2018\.  Graph convolutional neural networks for web-scale recommender systems. In _Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_.


# Observability of Forming Planets and their Circumplanetary Disks III. –
Polarized Scattered Light in Near-IR

J. Szulágyi1,2 & A. Garufi3,
1 Institute for Particle Physics and Astrophysics, ETH Zürich, Switzerland
2 Center for Theoretical Astrophysics and Cosmology, Institute for
Computational Science, University of Zürich, Switzerland
3 INAF, Osservatorio Astrofisico di Arcetri, Firenze, Italy
E-mail<EMAIL_ADDRESS>
(Accepted XX. Received XX; in original form 2019 May 17)

###### Abstract

There are growing amount of very high-resolution polarized scattered light
images of circumstellar disks. Nascent giant planets planets are surrounded by
their own circumplanetary disks which may scatter and polarize both the
planetary and stellar light. Here we investigate whether we could detect
circumplanetary disks with the same technique and what can we learn from such
detections. Here we created scattered light mock observations at 1.245 microns
(J band) for instruments like SPHERE and GPI, for various planetary masses
(0.3, 1.0, 5.0, 10.0 $\rm{M_{Jup}}$), disk inclinations (90, 60, 30, 0
degrees) and planet position angles (0, 45, 90 degrees). We found that the
detection of a circumplanetary disk at 50AU from the star is significantly
favored if the planet is massive ($\geq 5\rm{M_{Jup}}$) and the system is
nearly face-on ($\leq 30^{\circ}$). In these cases the accretion shock front
on the surface of the circumplanetary disks are strong and bright enough to
help the visibility of this subdisk. Its detection is hindered by the
neighboring circumstellar disk that also provides a strong polarized flux.
However, the comparison between the $PI$ and the $Q_{\phi}$ maps is a viable
tool to pinpoint the presence of the circumplanetary disk within the
circumstellar disk, as the two disks are behaving differently on those images.

###### keywords:

planets and satellites : detection – hydrodynamics – radiative transfer –
techniques: polarimetric

††pagerange: Observability of Forming Planets and their Circumplanetary Disks
III. – Polarized Scattered Light in Near-IR–9††pubyear: 2019

##  1 Introduction

Young, forming giant planets are surrounded by their circumplanetary disks
(CPDs), where their regular satellites will form eventually. Regardless
whether the planet formed via core accretion or disk instability scenario, the
circumplanetary disk forms in the last phase of the planet formation (e.g.
Kley, 1999; Lubow et al., 1999; Shabram & Boley, 2013; Ayliffe & Bate, 2009a,
b).

The circumplanetary disks characteristics have been studied numerically since
two decades now (e.g. Kley, 1999; Lubow et al., 1999; D’Angelo et al., 2003;
Ayliffe & Bate, 2009a, b; Machida et al., 2010; Ayliffe & Bate, 2012; Tanigawa
et al., 2012; Shabram & Boley, 2013; Gressel et al., 2013; Szulágyi et al.,
2014; Fujii et al., 2014; Perez et al., 2015; D’Angelo & Podolak, 2015;
Tanigawa et al., 2014; Zhu, Ju & Stone, 2016). As circumstellar disk, they
widely vary in mass from $10^{-4}$ $\rm{M_{planet}}$ (D’Angelo et al. 2003; G.
D’Angelo private communication; Gressel et al. 2013) till $\sim$
$\rm{M_{planet}}$ (Shabram & Boley, 2013; Szulágyi et al., 2017a) if the
planet formed via gravitational instability. Because the CPD is constantly fed
from the circumstellar disk (e.g. Szulágyi et al., 2014; Fung & Chiang, 2016),
the mass of the CPD depends on the circumstellar disk mass (apart from the
planetary mass Szulágyi 2017). The CPD’s temperature will depend on the semi-
major axis of the planet, the local density and opacity, the mass of the
planet and its age, the viscosity of the gas, among other factors (D’Angelo et
al., 2003; Papaloizou & Nelson, 2005; Ayliffe & Bate, 2009b, 2012; Gressel et
al., 2013; Szulágyi, 2017). The temperature of the CPD vary from thousands of
Kelvins in the forming planet vicinity to a few hundred Kelvins in the outer
subdisk. Of course, the CPD evolves in time, similarly to the circumstellar
disk, getting lighter and cooler during its lifetime (Szulágyi, 2017).

The characteristics of the circumplanetary disks are affecting their
detectability, hence creating mock observations from the hydrodynamical
simulations is a useful tool to plan and interpret real observations. Planet-
disk interactions, such as gaps has been studied on synthetic images (Dipierro
et al., 2015; Szulágyi et al., 2018; Zhang, et al., 2018). Circumplanetary
disks had been predicted to be seen with ALMA and VLA (Szulágyi et al., 2018;
Isella & Turner, 2018; Zhu, Andrews & Isella, 2018). Mock images of polarized
light about circumstellar disks helped us understanding what polarized light
observations can reveal about the circumstellar disk characteristics (Dong et
al., 2012). Synthetic observations of scattered light shed light on how
planet-disk interactions – especially spirals – are expected to look like
(Dong et al., 2015a, b; Fung & Dong, 2015; Dong, Fung, & Chiang, 2016). It has
also been suggested, that polarized light from the circumplanetary disk dust
could be detected in favorable circumstances (Stolker et al., 2017). The
detectability of the CPD in near-infrared as well as via spectral energy
distributions was investigated in Szulágyi et al. (2019). Actual detections of
circumplanetary disks are just began in 2019: two possible candidates from
ALMA continuum emission in the PDS70 disk (Isella et al., 2019), as well as
infrared excess detection of circumplanetary material around PDS70b
(Christiaens et al., 2019). With these few possible detections, the
characterization of circumplanetary disks from observations are still not yet
possible.

Unlike circumplanetary disks, circumstellar disks have been thoroughly
characterized from observations during the last decade thanks to optical/near-
IR instruments like VLT/SPHERE and GPI (e.g., Garufi et al., 2017; Rapson et
al., 2015) and to the (sub-)mm interferometer ALMA (e.g., Andrews et al.,
2018). Among the near-IR observations, the most successful technique to
directly image circumstellar disks is currently the polarized differential
imaging (PDI, Kuhn et al., 2001; Apai et al., 2004). This technique allows a
very good removal of the strong stellar flux by separating the polarized light
(mostly scattered light from the disk) from the unpolarized light (mainly
stellar light). Therefore, most of the available high-resolution near-IR maps
of circumstellar disks trace the polarized scattered light from the disk
surface. In principle, these polarized scattered light observations also open
the way to detect the circumplanetary disk the same way although this is yet
to be proven observationally.

In this paper we combine temperature-included (i.e. radiative) 3D gas
hydrodynamic simulations, with Monte-Carlo radiative transfer to create mock
observations about detecting the circumplanetary disk in scattered light with
and without polarization. In the first paper of this series, we looked at the
circumplanetary disk observability in sub-mm/radio wavelength (Szulágyi et
al., 2018). In the second paper, we reviewed the case for near-infrared and
spectral energy distributions (Szulágyi et al., 2019). In Szulágyi & Ercolano
(2020) we made predictions of hydrogen recombination line fluxes with
extinction and determined the planet-mass/planet accretion versus H-alpha,
Paschen-beta, Brackett-gamma line luminosity relationships.

##  2 Methods

We had a three step process for creating the mock images presented in this
work. First, we run 3D radiative hydrodynamic simulations of the circumstellar
disk with a forming planet embedded within (Sect. 2.1). Then we used the
RADMC-3D radiative transfer tool to create wavelength-dependent images of the
systems on 1.245 microns with polarization (Sect. 2.2). Finally, we convolved
the images with a diffraction limited PSF for the VLT/SPHERE instruments and
created polarization maps (Sect. 2.3).

###  2.1 Hydrodynamic Simulations

The hydrodynamic simulations in this study are the same as in our previous
paper (Szulágyi et al., 2019) of the series. In brief, we had a circumstellar
disk with a mass of $\sim 10^{-2}\mathrm{M_{Sun}}$ between 20 and 120 AU
around a solar-mass star, where a planet is forming at 50 AU. In four
different simulations, the planet masses were chosen to be a Saturn-mass, 1
Jupiter-mass, 5 Jupiter-masses and 10 Jupiter-masses (i.e. only one planet
present in each hydrodynamic run). We used the JUPITER code to carry our the
hydrodynamic calculations, that was developed by F. Masset and J. Szulágyi
(Szulágyi et al., 2014, 2016a) that not only solves Euler equations but also
the radiative transfer with the flux-limited diffusion approximation (two-
temperature approach Kley 1989, e.g.; Commerçon et al. 2011, e.g.). The
heating processes include adiabatic compression, viscous heating, shock
heating and stellar irradiation, while the cooling processes are adiabatic
expansion and radiative diffusion. The main source of heating in the
circumplanetary disk is the accretion process (Szulágyi et al., 2016a), as the
gas tries to fall onto the planet, leading to adiabatic compression in this
region that heats up the compressing gas. Furthermore, the accretion shock
front on the circumplanetary disk surface (Szulágyi & Mordasini, 2017) in case
of the higher mass planets are also strongly heated up. Viscous heating in the
CPD is secondary, with the chosen low viscosity value. Stellar irradiation is
computed based on solar flux, but its effect on the CPD is negligible. The
viscosity was a constant kinematic viscosity of
$10^{-5}\mathrm{a_{p}}^{2}\Omega_{p}$, where $\mathrm{a_{p}}$ is the semi-
major axis and $\Omega_{p}$ denotes the orbital frequency of the planet.

Given that we were particularly interested in the circumplanetary region,
where high-resolution is necessary to get the disk characteristics (density,
temperature, velocities) right, we used mesh refinement in this region. This
meant that while the circumstellar disk has been simulated with a lower
resolution (680 cells azimuthally over $2\pi$, 215 cells radially between 20
and 120 AU and 20 cells in the co-latitude direction over 7.4 degrees opening
angle from the midplane), the Hill-sphere of the planet were well resolved
with four levels of refinement. Each level doubled the resolution in each
spatial direction, hence the final resolution in the planet vicinity was 0.029
AU.

While the dust was not explicitly simulated within the hydrodynamics, its
effect on the temperature of the disk is taken into account through the dust
opacities (with the limit of assuming a constant dust-to-gas ratio of 1%). The
opacity table was equivalent to what was used in Szulágyi et al. (2019), and
included both gas and dust opacities.

###  2.2 RADMC-3D post-processing

RADMC-3D (Dullemond, 2012)111http://www.ita.uni-
heidelberg.de/~dullemond/software/radmc-3d/ radiative transfer tool was used
to create wavelength-dependent intensity images from the hydrodynamic
simulations. We used $5\times 10^{7}$ photons for these Monte-Carlo runs. We
run the RADMC-3D with the flux conservation option, which makes sure that the
total flux of the images are conserved as well, regardless the image
resolution. We ran two sets of images:

  * •

5000x5000 pixel resolution image on the entire circumstellar disk; ran
multiple times with planet positions of 0 deg, 45 deg, 90 deg, and
inclinations of 0 deg, 30 deg, 60 deg and 90 deg.

  * •

1000x1000 pixel images on the Hill-sphere (using zoomau command); ran 20 times
with randomly changing seed number and averaged at the end. We calculated the
variance between these 20 runs and verified that it is near zero, so
convergence was reached. These CPD region images were produced for the four
different inclinations (0 deg, 30 deg, 60 deg and 90 deg).

The dust-density files were created from the gas density (i.e. assuming that
these micron-sized dust grains are strongly coupled to the gas), by
multiplying the gas density in each cell with the dust-to-gas ratio of 1%. We
assumed thermal equilibrium, hence we used the dust temperature to be equal to
the gas temperature, except that the dust evaporation above 1500 K was taken
into account. This meant that in the cells hotter than this limit, the dust
density was set to zero, to be consistent with the radiative hydrodynamic
simulation, where the opacity table contains the dust evaporation as well. We
used the hydrodynamic simulation calculated temperature for the dust, because
that includes shock-heating, accretional heating, viscous heating which are
very important in the circumplanetary disk region and result in a hot planet
vicinity. We compared our results based on another method, where we used
thermal Monte Carlo simulation to calculate the dust temperature, using
radmc3dmctherm, but this method does not account for the main heating
mechanisms that take place in the planet vicinity, hence result in different
mock images (see discussion in Sect. 4).

The distance of the circumstellar disk was assumed to be 100 parsec for all
the calculations.

The hydrodynamic simulations cannot handle well optically thin, low-density
regions of the circumstellar disk, such as the disk atmosphere. In the hydro
simulations the disk opening angle was only 7.4 degrees, but real
circumstellar disks have a larger opening angle. Therefore we had to extend
the circumstellar disk in the vertical direction using an extrapolation
technique, before we ran the RADMC-3D calculations. The extrapolation was as
follows. First, we fitted Gaussian-functions to the density field in each cell
column (z-direction) separately, so that the vertical extent of the
circumstellar disk was 2.5 times larger than the original hydro simulation’s.
Second, in this circumstellar disk atmosphere region, we kept the temperature
as it is in the last (optically thin) co-latitude cell. Here the temperature
is high due to stellar irradiation, much higher than in the bulk of the
circumstellar disk (midplane regions). This meant that the temperature in the
circumstellar disk atmosphere was constant with co-latitude.

The dust opacities were identical to what had been used in Pohl et al. (2017).
It was assumed to be a mixture made of silicates (Draine 2003), carbon (Zubko
et al. 1996), and water ice (Warren & Brandt 2008) with fractional abundances
of 7%, 21%, and 42%, consistent with Ricci et al. (2010). The remaining 30%
was vacuum. The opacity of the mixture was determined by means of the
Bruggeman mixing formula. The absorption and scattering opacities,
$\kappa_{\mathrm{scat}}$ and $\kappa_{\mathrm{abs}}$, as well as the
scattering matrix elements $Z_{ij}$ were calculated for spherical, compact
dust grains with Mie theory considering the BHMIE code of Bohren & Huffman
(1984). The grain sizes were between 0.01 micron and 150 micron, with a power-
law index of -3.5.

###  2.3 Polarization maps

To compare our simulations to the available observations, we first convolved
the images with a rotationally symmetric 2D Gaussian Point-Spread-Function,
with a Full-width-half-maximum to be $1.22\cdot\lambda/D$, where $\lambda$ is
the wavelength and $D$ is the mirror-size of 8.2 meters (equivalent of VLT
mirror diameter). The RADMC-3D provides the set of Stokes parameters
$I,Q,U,V$. The polarized intensity map $PI$ was obtained through:

$PI=\sqrt{Q^{2}+U^{2}}$ (1)

An alternative treatment of the Stokes parameters is commonly used in
observational work, that is the creation of the tangential (sometimes called
radial or polar) parameters $Q_{\phi}$ and $U_{\phi}$ (Canovas et al., 2015;
Monnier et al., 2019). These are defined as:

$\begin{split}Q_{\phi}&=+Q\cos\,2\phi+U\sin\,2\phi\,,\\\
U_{\phi}&=-Q\sin\,2\phi+U\cos\,2\phi\,\end{split}$ (2)

with $\phi$ being the angle with respect to the stellar position (x0,y0)
calculated as:

$\phi=\mathrm{arctan}\frac{x-x_{0}}{y-y_{0}}\,$ (3)

By construction, $Q_{\phi}$ corresponds to $PI$ in the scenario of perfectly
centro-symmetric scattering and single scattering, whereas $U_{\phi}$ is
ideally expected to only contain noise.

##  3 Results

The obtained $I,\ PI,\ Q_{\phi}$ and $U_{\phi}$ maps of on the different
inclinations (0, 30, 60, 90 degrees), and different planetary positions (0,
45, 90 degrees). The J-band images of zero inclination are shown in Fig. 1,
while the other inclinations are in the Appendix A. Fig. 1 compares the
simulations of the four planetary masses considered: 0.3 $\rm{M_{Jup}}$, 1
$\rm{M_{Jup}}$, 5 $\rm{M_{Jup}}$, 10 $\rm{M_{Jup}}$. On the images, the planet
(and circumplanetary disk) always lies to the East at 50 AU from the central
star.

Figure 1: Polarized scattered light images at 1.245 microns (J band) with 0∘
inclination for the 10, 5, 1, and 0.3 Jupiter-mass cases (from top to bottom).
The columns are the $I$, $PI$, $Q_{\phi}$, $U_{\phi}$, and a zoom of the
$Q_{\phi}$ (top) and $U_{\phi}$ maps (bottom) on the CPD region, respectively.
The $PI$, $Q_{\phi}$, and $U_{\phi}$ images have the same color stretch to
highlight their relative flux brightness, whereas the zoomed maps have a
harder stretch with negative values shown in black. The white box indicates
the zoomed area of the last column. The yellow and red lines on the last
$Q_{\phi}$ map highlight the region used to calculate the contrast of CSD and
CPD, respectively. The assumed distance is 100 pc.

From these images, the main circumstellar disk (CSD) is always very bright in
$PI$ and its morphology resembles that of the $I$ images. The circumplanetary
disk is visible in the first two cases only, that is with a planet of 10 and 5
$\rm{M_{Jup}}$. Similar considerations apply to the $Q_{\phi}$ images and
these maps look very similar to the $PI$. On the other hand, the $U_{\phi}$
images do not show any significant signal except around the circumplanetary
disk in the first case. The reason for the CPD visibility only in the high-
mass planet cases, because in these cases the accretion shock front on the
surface of the CPD is strong (the velocity of the incoming accretion flow is
super-sonic), creating a hot, bright surface (Szulágyi & Mordasini, 2017). In
the smaller mass planet cases, the accretion flow is slower and sub-sonic,
hence the shock is not as strong and not as hot (see 4 for further discussion
on this point). This accretion shock is created by the incoming circumstellar
disk material, through the meridional circulation, a mass transfer between the
circumstellar and circumplanetary disks (Szulágyi et al., 2014; Fung & Chiang,
2016). The accretion shock on the circumplanetary disk surface hence
contributes to the observability and observational appearance of the CPD.

What is described above for the 0 inclination case (Fig. 1) also applies to
the other images created for the other inclinations (see Appendix A). The only
obvious differences are that the circumplanetary disk becomes decreasingly
evident with increasing inclination, and that some signal is recovered from
the $U_{\phi}$ image when the inclination is high, in agreement with the
theoretical prediction by Canovas et al. (2015).

Figure 2: Radial profile of the polarized-to-stellar light contrast of $PI$,
$Q_{\phi}$ and $U_{\phi}$ as described in Sect. 3.1. The profiles obtained
from three different planetary masses are displaced along the y-axis for a
better visualization. The vertical line indicates the planet location. The
$PI$ and $Q_{\phi}$ profiles are coincident except where the 5 $M_{\rm jup}$
mass planet is detected (see Sect. 3.2).

###  3.1 Polarized contrast

In this section, we provide a more quantitative analysis of the maps in Fig.
1, as well as of those shown in Appendix A. Measuring the amount of scattered
light from real observations is a challenging task because of the difficulties
in flux-calibrating the images and because the disk flux is directly dependent
on the stellar flux. Some authors (e.g., Avenhaus et al. 2018; Garufi et al.
2018) quantified the near-IR polarized light from the disk in relation to the
stellar flux, thus as to alleviate the dependence on the stellar brightness.
In particular, a way to do it is by dividing the observed polarized flux at a
certain disk radius, $F_{\rm pol}(r)$, by the stellar flux incident on that
disk region, $F_{*}/4\pi r^{2}$. This number contains information on both the
intrinsic albedo of particles (see e.g., Mulders et al. 2013) and on the
fraction of photons scattered toward the observer (see e.g., Stolker et al.
2016) and is thus sometimes referred to as (polarized) geometric albedo, or
contrast. This measurement is available for a relatively large number of real
circumstellar disks (see Garufi et al., 2017, 2018). From our simulations, we
obtained the aforementioned contrast along a 4 au-large radial cut oriented
toward the planet location. This profile is obtained from the $PI,\ Q_{\phi}$
and $U_{\phi}$ images as well, and is shown in Fig. 2 for some illustrative
cases. $Q_{\phi}$ almost always lies upon the profile of total polarized
intensity, and at the planet location there is zero excess from the presence
of the CPD in the case of Jupiter- and Saturn mass planets. For the higher
mass planets ($\geq$ 5 $\rm{M_{Jup}}$), at the CPD location there is a
contrast of 0.15 and 0.6 for $Q_{\phi}$ and $PI$, respectively.

A locally different flux recorded in the $Q_{\phi}$ and $PI$ maps indicates
that the pattern of the polarization diverges from centro-symmetric in the
circumstellar disk (see also Sect. 3.2). This can be appreciated by plotting
the polarization angles $\psi=0.5*\arctan(U/Q$) on top of the $Q_{\phi}$ maps,
as done in Fig. 3, zoomed to the Hill-sphere region. From the image, it is
clear how the CPD changes the polarization vectors locally in the 5 and 10
Mjup planet cases. Furthermore, in the 10 Mjup case the deviation include not
only the CPD ($\sim$ half of Hill-sphere), but the surrounding area as well
(the spiral wakes).

Figure 3: Polarization vectors overplotted on the $Q_{\phi}$ maps of the close
region around the CPD. In the first two cases, a centro-symmetric pattern is
visible. In the last two cases, polarization vectors are vertical, in line
with the local pattern of the protoplanetary disk.

We also extracted the contrast from the circumstellar and from the
circumplanetary disk by averaging the contribution from their respective
regions (0.04” for the CPD, 0.25” for the CSD that is exactly within the inner
and outer edge of the disk). The values thus obtained for the circumstellar
disk from the different simulations are comprised in a narrow interval of
values (from 1.5% to 3.2%). Compared to real disks, these numbers are
realistically high since the brightest disks ever observed in PDI have it up
to $\sim 2\%$ (see Garufi et al. 2017). This shows that if the CPDs have a
strong enough accretion-shock surface on them, they could be surprisingly
bright in some cases.

On the other hand, the contrast obtained around the circumplanetary disk span
enormously (from 4500% to $\lesssim 0.1\%$). From the $PI$ image, the contrast
of the 10 $\rm{M_{Jup}}$ case is always larger than 1 (i.e., more photons than
those incident from the star are detected) indicating a strong additional
source of photons to be scattered (the hot circumplanetary disk shock
surface). This observational scenario would be by itself a natural, robust
evidence of circumplanetary disk. However, for all the other planetary mass
cases that we studied the detection of the circumplanetary disk in polarized
light is less straightforward. Observationally, we can define a formal
threshold of 0.1% below which the signal is mostly noise (Garufi et al.,
2017). According to this criterion, 7 of the remaining 9 cases (3 planet
masses, 3 inclinations) should still be regarded as detection. We caution,
however, that our images does not contain extra source of noise, that could
make the detections even more difficult as described here. So these detection
numbers should be regarded as idealistic, best-case scenarios.

We must nonetheless consider the effect of the circumstellar disk itself that
may still be present at the planet location (in particular for the 0.3
$\rm{M_{Jup}}$ case where the disk gap is more shallow than for the more
massive planets) and leaves the same imprint on the scattered-light images. In
this regard, we noticed that the contrast around the planet decreases toward
smaller masses but then increases again for the lowest-mass case.

###  3.2 Circumstellar versus circumplanetary disk signal

Our simulations show that it is formally possible to distinguish between the
scattered light from the circumstellar and from the circumplanetary disk by
comparing the contrast from the $PI$ and $Q_{\phi}$ images. In fact, for the
two largest-mass planet scenarios the polarized contrast around the planet
calculated from these two maps significantly differ (up to a factor 50)
whereas in the 1.0 $\rm{M_{Jup}}$ case only a minor ratio ($\sim 30\%$) is
visible, and no difference is appreciable in the 0.3 $\rm{M_{Jup}}$ case.
Conversely, for all our simulations the circumstellar disk signal from the
$PI$ and $Q_{\phi}$ images is very similar (always within $10\%$). This
behaviour can be appreciated from Fig. 4. Strong discrepancies between $PI$
and $Q_{\phi}$ are expected when the scattered light deviates from a centro-
symmetric pattern (see Fig. 3), which is the assumption under which $Q_{\phi}$
is constructed (see Eq. 2). In the presence of a circumplanetary disk, photons
are not expected to be scattered in such a pattern since the star is no longer
the only source of photons (see on Fig. 1 for 10 Jupiter-mass case).
Therefore, the comparison of the polarized contrast from the $PI$ and
$Q_{\phi}$ images is a simple but potentially powerful manner to discriminate
the presence of a circumplanetary disk.

Figure 4:  $PI$/$Q_{\phi}$ contrast ratios of circumstellar- and
circumplanetary disks for different planet masses and different inclinations.
From these simulations, in the 10 and 5 $\rm{M_{Jup}}$ cases, as well as
marginally in the 1 $\rm{M_{Jup}}$ case, the signal from the circumplanetary
disk can be observationally disentangled from the circumstellar disk signal.

Comparing the $PI$/$Q_{\phi}$ contrast ratios of the CPD with various position
angles of the planet (0, 45, 90 degrees) does not show a clear trend towards
either direction, hence we only show the position angle 0 on Fig. 4. However,
for all position angles considered, we found that the $PI$/$Q_{\phi}$ contrast
ratio decreases with increasing inclination angle.

##  4 Discussion

The results of Sect. 3 suggest that the parallel employment of $PI$ and
$Q_{\phi}$ maps and of multiple wavebands may help reveal a circumplanetary
disk. In fact, while the light scattered off by the CPD could easily be
confused for a substructure of the circumstellar disk, its polarization
pattern is different. This results in a local divergence between the
tangential and the total component of the polarized light, as indicated by the
$Q_{\phi}$ and $PI$ maps respectively. In this work, we showed that this
effect is appreciable for massive planets, with mass of the of at least
$5\,M_{\rm jup}$.

As of today, the paucity of planets detected in circumstellar disks does not
allow to test these predictions. Nonetheless, the ideal threshold of
$5\,M_{\rm jup}$ over which the scattered light from a CPD becomes detectable
is close to the common mass upper limit determined in a number of objects
(Claudi et al., 2019; Maire et al., 2017; Mesa et al., 2019a, b) indicating
that this approach could in principle be used in parallel to the typical
differential techniques to detect the planetary thermal light.

The results presented here are greatly affected by the temperature
calculation. We used the hydrodynamic simulation calculated temperatures (via
flux limited diffusion approximation), where assumed perfect thermal
equilibrium between the micron sized dust and the gas (Fig. 5 top-left panel).
The more traditional method to calculate the dust temperature is thermal Monte
Carlo computation, that can be done with RADMC-3D’s mctherm command (Fig. 5
top-right panel; calculated with $10^{6}$ photon packages). However, latter
method does not include the accretional heating due to adiabatic compression,
shock-heating due to accretion shock fronts, and viscous heating, all which
heats up the planet vicinity. Hence, the circumplanetary region does not show
up on the mctherm temperature maps (Fig. 5 top-right panel), which also means
that this region would not be very visible on mock observations. The
difference of the two temperature calculation are shown on Fig. 5 bottom
panels. Everywhere else in the simulation box (i.e. in the circumstellar disk)
the temperature difference between the two method is always smaller than
10-20K, the only region where the difference is at minimum 50 K (but ranges up
to 1600 K difference) is the accreting planet vicinity. This test shows that
for creating mock observations, the way the temperature is calculated is very
important for the outcome.

Figure 5:  Comparison of the temperature calculations in the 10 Jupiter-mass
case. Top-left panel: hydrodynamic simulation calculated temperature in the
midplane that include adiabatic compression, shock heating and viscous
heating. Top right panel: RADMC-3D calculated temperature field in midplane
with thermal Monte-Carlo (mctherm). Bottom-left panel: difference of the
temperature between the hydrodynamic simulation calculated temperature and
RADMC-3D’s Monte-Carlo calculated one, clearly there is a large difference in
the circumplanetary region. Bottom-right panel: the region in the entire
simulation box, where the temperature difference between the two methods are
at least 50 K – clearly this is only the circumplanetary region (a vertical
slice is shown of this area). The planet and the inner circumplanetary disk is
the orange region in the middle, the shock front on the circumplanetary disk
created by the accretion stream is also visible as a horizontal yellow line
above the planet.

In the 5 and 10 Jupiter-mass planet cases the circumplanetary region pops up,
and some spurious photons can be seen to originate from this region (Fig. 1).
Closer inspection showed that the accretion shock front created on the surface
of the circumplanetary disk (Fig. 6) due to the incoming mass influx from the
circumstellar disk (via the meridional circulation; Szulágyi et al. 2014) is
the origin of the noisiness on the mock images on Fig. 1. This shock front is
strong and hot ($>1000$ K) enough only in the 5 and 10 Jupiter-mass cases, not
for the smaller planets. This is due to the the fact that the incoming
accretion flow velocity scales with the planetary mass (nearly free-fall
velocity). The higher velocity of the influx creates stronger and hotter shock
fronts on the circumplanetary disk surface, with increasing planetary mass
(Szulágyi, 2017). This shock-front is helping the observability of the
circumplanetary disk, hence it is also important to self-consistently include
it when creating mock observations on circumplanetary disks.

Figure 6:  Hydrodynamic simulation calculated temperature field (threshold) in
the circumplanetary disk (10 Jupiter-mass case), where the accretion shock-
front ($>1000$ K) on the disk surface is visible by yellow, orange and red
areas above the planet. The noisiness of the circumplanetary region in the
mock observations (e.g. Fig. 1) is due to this shock front. However, this
shock front is also helping the detection of the circumplanetary disk.

The results also depend on the optical depth. In this work we took care of
dust evaporation above silicate evaporation temperature (1500 K), which meant
that the dust density were put to zero, where the temperatures were rising
above this limit (this region is near the planet, e.g. in the shock front on
the circumplanetary disk surface). We made a test run for the 10 Jupiter-mass
planet case (Fig. 7), where we compared cases with and without dust
evaporation taken into account, for the Monte-Carlo temperature calculation
versus the hydrodynamic simulation resulted temperatures. Fig. 7 shows that
when the dust evaporation is not included, the circumplanetary disk region
(and its surface shock front) does not show up on the mock images. This is
partially due to the different optical thickness – if evaporation taken into
account, and dust is no longer contribute to the local optical depth, the
region became more optically thin (gas-only case). This again highlights that
the modeling of the planetary vicinity should be carefully done to incorporate
various physical effects when making observational predictions.

Figure 7:  Comparison between the RADMC-3D total intensity images when the
dust temperature was calculated with the hydrodynamic simulation radiative
transfer module (that includes shock heating, and accretional heating due to
adiabatic compression, and viscous heating) - top row, and with RADMC-3D’s
mctherm Monte-Carlo method - bottom row. In the left column the dust
evaporation was not taken into account, while on the right evaporation was
included (i.e. the dust density was set to zero where the temperature rose
above 1500 K). This test run was done with $10^{6}$ photon packages on 1.245
microns, on the 10 Jupiter-mass simulation, units are in Jy/pixel, scaled to
100 pc distance.

Our models only covers part of the parameter space. We assumed a fixed dust-
to-gas ratio of 0.01 even in the circumplanetary disk (Dra̧żkowska & Szulágyi,
2018), however real disks can have smaller and larger values than this (e.g.
Youdin & Goodman, 2005; van der Marel et al., 2013; Dra̧żkowska & Dullemond,
2014; Birnstiel et al., 2012; Williams & Best, 2014; Ansdell et al., 2016),
which might affect the results. In this work we have considered the planets to
be 50 AU from the star, but the circumplanetary disk-circumstellar disk
contrast can be very different if the circumplanetary disk at another
distance. Circumplanetary disks closer to the star tend to be more optically
thick, and hotter than the more distant ones.

For the circumstellar disk mass we considered an average value of 0.01
$\rm{M_{Sun}}$, and the radial extent was between 20-120 AU, similar to a
transitional disk with an inner cavity. While the circumplanetary disk mass
linearly scales with the circumstellar disk mass (Szulágyi, 2017), the changes
in mass will also result in different optical depth, which can affect the
results described here. The large, optically thin inner 20 AU can also affect
the results.

The hydrodynamic simulations did not include magnetic fields, e.g. the fields
of the disks, which might affect the dust density distribution (Gressel et
al., 2013). Since the planet interior structure was not part of the
simulation, any temperature in the planet region is a lower limit to the
radiation from the planet. This will affect the radiation field in/from the
CPD, as well as its spatial structure.

##  5 Summary

In this work we investigated polarized scattered light detectability of
circumplanetary disks surrounding nascent planets. We ran hydrodynamic
simulations with mesh refinement to resolve sufficiently the circumplanetary
disk. We used radiative transfer included hydrodynamics, to realistically
estimate the temperature. Then, we post-processed the simulations with
RADMC-3D Monte-Carlo radiative transfer software to create polarized light
images in J band. We added convolution with a PSF-size at the diffraction
limit, assuming 8.2 meter mirror, like VLT.

We considered different planetary mass cases: Saturn, 1, 5, 10 $\rm{M_{Jup}}$,
different disk inclinations of 0, 30, 60, 90 degrees, and various planetary
locations (0, 45, 90 degrees). The planets were embedded in a 0.01
$\rm{M_{Sun}}$ circumstellar disk, 50 AU away from their star, which was
assumed to be a Sun-equivalent.

Our $I,\ PI,\ Q_{\phi}$ and $U_{\phi}$ images revealed that the
circumplanetary disk detection is only possible in the case of very massive
planets (5 and 10 $\rm{M_{Jup}}$), although it is highly dependent on how
optically thick is the circumplanetary disk (i.e. how much dust it contains,
and what is the temperature there, whether the dust is evaporated or not). In
these high mass planet cases the accretion shock front on the surface of the
circumplanetary disk (Szulágyi & Mordasini, 2017) is so strong and luminous,
that it helps the observability of this subdisk. Therefore the inclusion of
this shock-front is important when modelling observability of the
circumplanetary disk.

The circumplanetary disk detection is challenging in polarized light, not only
because of sensitivity but also due to the contrast with the circumstellar
disk. However, we showed that, ideally speaking, it is possible to distinguish
between the two disk’s contributions by comparing the total polarized light
(from the $PI$ image) and the centro-symmetric polarized light (from the
$Q_{\phi}$ image), as well as by finding stronger polarized colors in the
circumplanetary disk than in the neighboring circumstellar disk.

In conclusion, while circumplanetary disk detection might be challenging in
polarized light, the $PI$/$Q_{\phi}$ images can be possible tools to detect
the circumplanetary disk within the circumstellar disk.

## Acknowledgments

We thank for Adriana Pohl providing the opacity table, including the
polarization matrix. J.Sz. thanks for the financial support through the Swiss
National Science Foundation (SNSF) Ambizione grant PZ00P2_174115. Furthermore,
these results are part of a project that has received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (Grant agreement No. 948467). We also
acknowledge support from the project PRIN-INAF 2016 The Cradle of Life -
GENESIS-SKA (General Conditions in Early Planetary Systems for the rise of
life with SKA) and from INAF/Frontiera (Fostering high ResolutiON Technology
and Innovation for Exoplanets and Research in Astrophysics) through the
”Progetti Premiali” funding scheme of the Italian Ministry of Education,
University, and Research. Computations partially have been done on the ”Piz
Daint” machine hosted at the Swiss National Computational Centre and partially
carried out on ETH Zürich’s Euler machine.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Avenhaus et al. (2018) Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44
  * Ayliffe & Bate (2009a) Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 393, 49
  * Ayliffe & Bate (2009b) Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 397, 657
  * Ayliffe & Bate (2012) Ayliffe B. A., Bate M. R., 2012, MNRAS, 427, 2597
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41.
  * Ansdell et al. (2016) Ansdell, M., Williams, J. P., van der Marel, N., et al. 2016, ApJ, 828, 46
  * Apai et al. (2004) Apai, D., Pascucci, I., Brandner, W., et al. 2004, A&A, 415, 671
  * Birnstiel et al. (2012) Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148
  * Bohren & Huffman (1984) Bohren, C. F., & Huffman, D. R. 1984, Nature, 307, 575
  * Canovas et al. (2015) Canovas H., Ménard F., de Boer J., Pinte C., Avenhaus H., Schreiber M. R., 2015, A&A, 582, L7
  * Christiaens et al. (2019) Christiaens V., Cantalloube F., Casassus S., Price D. J., Absil O., Pinte C., Girard J., et al., 2019, ApJL, 877, L33. doi:10.3847/2041-8213/ab212b
  * Claudi et al. (2019) Claudi R., Maire A.-L., Mesa D., Cheetham A., Fontanive C., Gratton R., Zurlo A., et al., 2019, A&A, 622, A96. doi:10.1051/0004-6361/201833990
  * Commerçon et al. (2011) Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35
  * D’Angelo et al. (2003) D’Angelo, G., Henning, T., & Kley, W. 2003, ApJ, 599, 548
  * D’Angelo & Podolak (2015) D’Angelo, G., & Podolak, M. 2015, ApJ, 806, 203
  * Dipierro et al. (2015) Dipierro G., Pinilla P., Lodato G., Testi L., 2015, MNRAS, 451, 974
  * Dong et al. (2012) Dong, R., Rafikov, R., Zhu, Z., et al. 2012, ApJ, 750, 161
  * Dong et al. (2015a) Dong R., Zhu Z., Rafikov R. R., Stone J. M., 2015, ApJ, 809, L5
  * Dong et al. (2015b) Dong R., Hall C., Rice K., Chiang E., 2015, ApJ, 812, L32
  * Dong, Fung, & Chiang (2016) Dong R., Fung J., Chiang E., 2016, ApJ, 826, 75
  * Draine (2003) Draine, B. T. 2003, ApJ, 598, 1026
  * Dra̧żkowska & Szulágyi (2018) Dra̧żkowska J. & Szulágyi J., 2018, ApJ, 866, 142
  * Dra̧żkowska & Dullemond (2014) Dra̧żkowska, J., & Dullemond, C. P. 2014, A&A, 572, A78
  * Dullemond (2012) Dullemond, C. P. 2012, Astrophysics Source Code Library, ascl:1202.015
  * Fujii et al. (2014) Fujii, Y. I., Okuzumi, S., Tanigawa, T., & Inutsuka, S.-i. 2014, ApJ, 785, 101
  * Fung & Chiang (2016) Fung, J., & Chiang, E. 2016, arXiv:1606.02299
  * Fung & Dong (2015) Fung J., Dong R., 2015, ApJ, 815, L21
  * Garufi et al. (2017) Garufi A., et al., 2017, A&A, 603, A21
  * Garufi et al. (2017) Garufi, A., Benisty, M., Stolker, T., et al. 2017, The Messenger, 169, 32.
  * Garufi et al. (2018) Garufi A., et al., 2018, A&A, 620, A94
  * Gressel et al. (2013) Gressel, O., Nelson, R. P., Turner, N. J., & Ziegler, U. 2013, ApJ, 779, 59
  * Isella & Turner (2018) Isella A., Turner N. J., 2018, ApJ, 860, 27
  * Isella et al. (2019) Isella A., Benisty M., Teague R., Bae J., Keppler M., Facchini S., Pérez L., 2019, ApJL, 879, L25
  * Keppler, et al. (2018) Keppler M., et al., 2018, A&A, 617, A44
  * Kley (1989) Kley W., 1989, A&A, 208, 98
  * Kley (1999) Kley W., 1999, MNRAS, 303, 696
  * Kraus & Ireland (2012) Kraus, A. L., & Ireland, M. J. 2012, ApJ, 745, 5
  * Kuhn et al. (2001) Kuhn, J. R., Potter, D., & Parise, B. 2001, ApJ, 553, L189.
  * Lubow et al. (1999) Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001
  * Machida et al. (2010) Machida, M. N., Kokubo, E., Inutsuka, S.-I., & Matsumoto, T., 2010, MNRAS, 405, 1227
  * Maire et al. (2017) Maire A.-L., Stolker T., Messina S., Müller A., Biller B. A., Currie T., Dominik C., et al., 2017, A&A, 601, A134. doi:10.1051/0004-6361/201629896
  * Mesa et al. (2019a) Mesa D., Langlois M., Garufi A., Gratton R., Desidera S., D’Orazi V., Flasseur O., et al., 2019, MNRAS, 488, 37. doi:10.1093/mnras/stz1662
  * Mesa et al. (2019b) Mesa D., Keppler M., Cantalloube F., Rodet L., Charnay B., Gratton R., Langlois M., et al., 2019, A&A, 632, A25. doi:10.1051/0004-6361/201936764
  * Monnier et al. (2019) Monnier J. D., Harries T. J., Bae J., Setterholm B. R., Laws A., Aarnio A., Adams F. C., et al., 2019, ApJ, 872, 122
  * Mulders et al. (2013) Mulders G. D., Min M., Dominik C., Debes J. H., Schneider G., 2013, A&A, 549, A112
  * Müller et al. (2018) Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2
  * Murakawa (2010) Murakawa, K. 2010, A&A, 518, A63.
  * Papaloizou & Nelson (2005) Papaloizou J. C. B., Nelson R. P., 2005, A&A, 433, 247
  * Perez et al. (2015) Perez, S., Dunhill, A., Casassus, S., et al. 2015, ApJL, 811, L5
  * Pérez, et al. (2019) Pérez S., et al., 2019, arXiv, arXiv:1906.06305
  * Pohl et al. (2017) Pohl A., et al., 2017, A&A, 605, A34
  * Rapson et al. (2015) Rapson, V. A., Kastner, J. H., Millar-Blanchaer, M. A., et al. 2015, ApJ, 815, L26.
  * Ricci et al. (2010) Ricci L., Testi L., Natta A., Neri R., Cabrit S., Herczeg G. J., 2010, A&A, 512, A15
  * Sallum et al. (2015) Sallum, S., Follette, K. B., Eisner, J. A., et al. 2015, Nature, 527, 342
  * Schmid et al. (2006) Schmid H. M., Joos F., Tschan D., 2006, A&A, 452, 657
  * Shabram & Boley (2013) Shabram, M., & Boley, A. C. 2013, ApJ, 767, 63
  * Stolker et al. (2017) Stolker T., Min M., Stam D. M., Mollière P., Dominik C., Waters L. B. F. M., 2017, A&A, 607, A42
  * Stolker et al. (2016) Stolker T., Dominik C., Min M., Garufi A., Mulders G. D., Avenhaus H., 2016, A&A, 596, A70
  * Szulágyi et al. (2014) Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65
  * Szulágyi et al. (2016a) Szulágyi, J., Masset, F., Lega, E., et al. 2016, MNRAS, 460, 2853
  * Szulágyi (2017) Szulágyi, J. 2017, ApJ, 842, 103
  * Szulágyi et al. (2017a) Szulágyi, J., Mayer, L., & Quinn, T. 2017, MNRAS, 464, 3158
  * Szulágyi et al. (2018) Szulágyi, J., Plas, G. v. d., Meyer, M. R., et al. 2018, MNRAS, 473, 3573
  * Szulágyi et al. (2019) Szulágyi J., Dullemond C. P., Pohl A., Quanz S. P., 2019, MNRAS, 1269
  * Szulágyi & Mordasini (2017) Szulágyi, J. & Mordasini, C., 2017, MNRAS, 465, L64
  * Szulágyi & Ercolano (2020) Szulágyi J., Ercolano B., 2020, arXiv, arXiv:2002.09918
  * Tanigawa et al. (2012) Tanigawa, T., Ohtsuki, K., & Machida, M. N., 2012, ApJ, 747, 47
  * Tanigawa et al. (2014) Tanigawa, T., Maruta, A., & Machida, M. N. 2014, ApJ, 784, 109
  * van der Marel et al. (2013) van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199
  * Wagner et al. (2018) Wagner, K., Follete, K. B., Close, L. M., et al. 2018, ApJL, 863, L8
  * Warren & Brandt (2008) Warren, S. G., & Brandt, R. E. 2008, J. Geophys. Res., 113, D14220
  * Williams & Best (2014) Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59
  * Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459
  * Zhang, et al. (2018) Zhang S., et al., 2018, ApJ, 869, L47
  * Zhu, Ju & Stone (2016) Zhu Z., Ju W., Stone J. M., 2016, ApJ, 832, 193
  * Zhu, Andrews & Isella (2018) Zhu Z., Andrews S. M., Isella A., 2018, MNRAS, 479, 1850
  * Zubko et al. (1996) Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321

##  Appendix A Maps

Similarly to Fig. 1, Figs. 8, 9 show the $I,\ PI,\ Q_{\phi}$ and $U_{\phi}$
maps for the four planetary masses explored. In particular, Fig. 8 shows 30
deg inclination and Fig. 9 represents 60 deg inclination.

Figure 8: Same as Fig. 1, but for 30∘ inclination. Figure 9: Same as Fig. 1,
but for 60∘ inclination.