RHADaMAnTe: An Astro Code to Estimate the Energy Spectral Distribution of a Curved Wall of a Gap Opened by a Forming Planet in a Protoplanetary Disk

1. Introduction

In a small fraction of young stellar objects (YSOs) surrounded by disks, observations have discovered a low excess radiation in the near–infrared but a high excess in longer wavelengths. This has been interpreted as an evidence that these disks called transitional disks (TDs) have central holes which have practically no dust [1]. More recently, disks showing a significant excess in the near–infrared have been discovered. Such an excess indicates the presence of an optically thick inner disk. This inner disk is separated from an outer disk which also has a high optical depth. In this way the spectral energy distribution (SED) suggests the incipient development of a gap between both disks, these disks are called pre–transitional disks (Pre–TDs) [2]. Several physical mechanisms have been suggested to explain the gaps or holes in protoplanetary disks. The one implemented in this work is driven by forming giant planets.

A key element that produces characteristic features in the SEDs of protoplanetary disks is the outer wall of the gap or hole. To simplify SED wall models, it is often assumed that the wall is vertical and frontally irradiated by the central star [2,3]. But this assumption is physically wrong [4]. For dust sublimation walls, located near from the star, it has been proposed that the wall is curved, where the dust grain growth and its fall into the mid–plane of the disk, and the gas density high–dependence on sublimation temperature are the physical mechanisms responsible for such curvature [5].

In order to create synthetic SEDs of protoplanetary disks with gaps or holes having inner curved walls, I have developed a computational code called rhadamante. This code is based on an older code which suggests that the inner vertical wall of the outer disk can explain the mid–infrared spectrum of the low–mass pre main–sequence star CoKu Tau/4 [6].

To test the code, I present a model of a truncated dusty disk –a disk with an inner hole– that accounts for the Spitzer Infrared Spectrograph observations of the low-mass pre main–sequence star LkCa 15. In this model the mid–infrared spectral energy distribution (between 10 and 25µm) arises from the inner curved wall of the gap in the disk.

1.1. Dust in Protoplanetary Disks

Dust is a pretty important component of protoplanetary disks surrounding young stars. The growth of dust grains from sizes of microns to centimeters or larger grains is the first step in planet formation.

The dust grains in protoplanetary disks follow a size distribution based on a single power law $f\left(a\right)\propto a^p$ [[7], MNR], where the maximum dust size, $a_{max}$, the minimum dust size, $a_{min}$, and the power law index, p, are different for each grain species.

The dust composition in protoplanetary disks has been extensively studied through mid-IR observations. This dusty mixture includes silicates (mainly), carbonaceous grains, poly-cyclic aromatic hydrocarbons, and sulfide-bearing grains [8].

Silicate grains are the best understood dust component. The most abundant crystalline silicates are olivine, which is magnesium-rich (Fo₉₀), and pyroxene. Olivine series range from forsterite, Mg₂SiO₄ (denoted as Fo₁₀₀), to the fayalite, Fe₂SiO₄ (denoted as Fo₀). While pyroxene series range from the enstatite, MgSiO₃ (denoted En₁₀₀), to the ferrosilite, FeSiO₃ (denoted En₀).

Modeling of the observed spectra expects amorphous silicate grains to exist in protoplanetary disks [see, e.g. [9]; [10]]. The composition of these grains are glass with embedded metals and sulfides, and series ranging from ferromagnesian silica to Fe–Mg-bearing aluminosilica. These grains are difficult to observe directly from infrared spectroscopy. Their spectral signature observed is a combination of grain composition, shape, size, and structure, making difficult to isolate the pure amorphous silicate signal.

Carbonaceous grains, including amorphous and graphite elemental carbon, are difficult to detect in the infrared. However, grain modeling suggests these grains are needed in order to explain the observed infrared spectra of protoplanetary disks [9].

Nano–diamonds, from sizes of <1 nm to ∼ 10 nm, are found in protoplanetary disks. Diamond emission coming from the inner region of the disk (i.e. < 15 AU) at 3.43 and 3.53 $\mu$m has been detected in disks [see, e.g. [11]].

The presence of poly-cyclic aromatic hydrocarbons (PAHs) has been detected in the surface layers of some protoplanetary disks [see, e.g. [12]]. Disks surrounding higher-mass stars, such as Herbig stars, show more PAHs emission [13] than disks surrounding lower-mass stars, such as T Tauri stars [14]. Protoplanetary disks with a flaring outer surface show significantly more PAH emission [13]. It follows that PAHs exist in all disks, but they can only be detected, as infrared emission, when ultraviolet radiation from the central star is able to excite them. The discovery of weak PAH features in T Tauri stars supports this idea [see, e.g. [14]].

Other dusty components in protoplanetary disks are iron-nickel sulfides grains (FeS, NiS) and water ice (H₂O). Sulfide emision around 23 $\mu$m has been detected in the emission spectra of protoplanetary disks [15]. While water ice emission has been identified at 3$\mu$m [16], 44$\mu$m [17], 60$\mu$m [18] and 62$\mu$m [19].

1.2. LkCa 15

LkCa 15 is a K5-type [20] T Tauri star located in the nearby ($145\pm 15\mathrm{pc}$) the Taurus-Auriga Star Forming Region [21]. The mass of the central star is $0.97\pm 0.03M_{\odot }$[20], it has an effective temperature of 4370K [22] and a radius of $1.6R_{\odot }$[23]. Three planet candidates have been detected: LkCa 15b (semi major axis $a=14.7\pm 2.1\mathrm{AU}$) [24], LkCa 15c ($a=18.6\pm 2.5\mathrm{AU}$) and LkCa 15d ($a=18.0{\pm }_{5.4}^{6.7}\mathrm{AU}$), with masses lower than 5–10$M_{♃}$, for the two first planets, and $\le 0.5M_{♃}$ for the third one [25].

Observations of the far-ultraviolet (1100–2200Å) radiation field and the near–to mid–IR (3-13.5$\mu$m) spectral energy distribution of LkCa 15, from the Space Telescope Imaging Spectrograph (STIS) indicate the existence of an inner disk gap of a few astronomical units [26].

LkCa 15 has an inner disk, a gap and an outer disk [27]. Using the Spitzer data, LkCa 15 has been classified as a pre-transitional disk [28], and it has been showed that the inner hole is not devoid of dust between 0.1 and 5 AU. The outer disk extends from 46 to 800 AU [29].

Recent observations from Gemini NIRI suggest that a single massive planet would be capable of opening a gap as large as the one observed in the LkCa 15 disk [22]. This assumption leads me to use a mass of ∼$10M_{♃}$ for the planet candidate in the current work.

Models of the LkCa 15 disk SED show that the inner edge of the disk has a radius of ∼58 AU [23]; this outer disk has a mass of ∼$50{\mathrm{M}}_{♃}$ [30].

2. Geometry of the Wall Projected on the Sky

To find the two–dimensional geometry of a wall gap, I implement the artemise code [4]. Which is a computational and geometrical code that analyses a tri–dimensional simulation of the disk–planet interaction by considering the wall is located at the points $(x,y,z)$ where the disk optical depth is ${\tau }_{wall}=\frac{2}{3}$. Simulations are done with fargo-3D code [31] under some specific parameters of the young stellar object to be studied.

2.1. Inclined Walls

Definition 1.

An inclined two–dimensional wall is a line segment with boundaries $(R_{\mathrm{wall}}^{\mathrm{up}},H_{\mathrm{wall}}^{\mathrm{up}})$ and $(R_{\mathrm{wall}}^{\mathrm{down}},H_{\mathrm{wall}}^{\mathrm{down}})$, with $R_{\mathrm{wall}}^{\mathrm{up}}>R_{\mathrm{wall}}^{\mathrm{down}}$ and $H_{\mathrm{wall}}^{\mathrm{down}}=-H_{\mathrm{wall}}^{\mathrm{up}}$, as seen in Figure 1a. By rotating this line segment around z–axis, a conic ring lying in the Euclidean space $(x,y,z)$, is obtained, as seen in Figure 1b. This ring is the tri–dimensional conic wall.

In the coordinate system $(x,y,z)$ the star is centered at the origin, here z–axis is the disk rotation axis, and the plane $(x,y)$ is the disk mid–plane. For simplicity, I also consider the cylindrical coordinate system $(R,\theta ,z)$, such that all points on the wall superior boundary have coordinates z, $x=R_{\mathrm{wall}}^{\mathrm{up}}cos\theta$ y $y=R_{\mathrm{wall}}^{\mathrm{up}}sin\theta$, whereas all points on the wall inferior boundary have coordinates z, $x=R_{\mathrm{wall}}^{\mathrm{down}}cos\theta$ y $y=R_{\mathrm{wall}}^{\mathrm{down}}sin\theta$. Since the protoplanetary disk is assumed to be projected on the plane of the sky $(X,Y)$, as seen in Figure 2, I consider a third coordinate system $(X,Y,Z)$ also centered at the star, where the Z–axis is the line of sight. When the disk is face–on the coordinate systems coincide. There exists a transformation between the three coordinate systems:

$$X=x,$$

(1a)

$$Y=ycos\left(i\right)-zsin\left(i\right),$$

(1b)

$$Z=zcos\left(i\right)-ysin\left(i\right),$$

(1c)

where i is the disk inclination angle, that is, the angle between the z–axis and the plane of the sky $(X,Y)$.

The amount of visible surface of the wall, projected on the plane of the sky, depends on the disk inclination angle, and there are two possibilities: (i) when the star is visible (corresponding to $\delta <1$, see Equations (4) and (9) for a definition of $\delta$), and (ii) when the star in invisible (corresponding to $\delta >1)$, as seen in Figure 2. A surface element of the visible area is $dA=dXdY=-R_{\mathrm{wall}}sin\theta d\theta dY$, with $R_{\mathrm{wall}}^{\mathrm{down}}\le R_{\mathrm{wall}}\le R_{\mathrm{wall}}^{\mathrm{up}}$.

Let $\mathcal{R}$ be the visible surface of the wall projected on the plane $(X,Y)$ for both cases, as seen in Figure 3. Then the boundary of this region is defined by two ellipses ${\epsilon }_{\mathrm{up}}$ and ${\epsilon }_{\mathrm{down}}$ (see Appendix A.1) given by the projections of the up and down edges of the tri–dimensional conic wall. The up ellipse ${\epsilon }_{\mathrm{up}}$ is defined as $Y_{\mathrm{up}}=\mathrm{s}{Y}_{\mathrm{up}}\cup \mathrm{i}{Y}_{\mathrm{up}}$, where $\mathrm{s}{Y}_{\mathrm{up}}$ and $\mathrm{i}{Y}_{\mathrm{up}}$ are the superior and inferior parts of this ellipse, respectively, such that

$${\displaystyle \frac{\mathrm{s}{Y}_{\mathrm{up}}}{R_{\mathrm{wall}}^{\mathrm{up}}}}=cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}^{\mathrm{up}}}}\right)}^2}+{\displaystyle \frac{H_{\mathrm{wall}}^{\mathrm{up}}}{R_{\mathrm{wall}}^{\mathrm{up}}}}sin\left(i\right),$$

(2a)

$${\displaystyle \frac{\mathrm{i}{Y}_{\mathrm{up}}}{R_{\mathrm{wall}}^{\mathrm{up}}}}=-cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}^{\mathrm{up}}}}\right)}^2}+{\displaystyle \frac{H_{\mathrm{wall}}^{\mathrm{up}}}{R_{\mathrm{wall}}^{\mathrm{up}}}}sin\left(i\right).$$

(2b)

Similarly, the down ellipse ${\epsilon }_{\mathrm{down}}$ is defined as $Y_{\mathrm{down}}=\mathrm{s}{Y}_{\mathrm{down}}\cup \mathrm{i}{Y}_{\mathrm{down}}$, where $\mathrm{s}{Y}_{\mathrm{down}}$ and $\mathrm{i}{Y}_{\mathrm{down}}$ are the superior and inferior parts of this ellipse, respectively, such that

$${\displaystyle \frac{\mathrm{s}{Y}_{\mathrm{down}}}{R_{\mathrm{wall}}^{\mathrm{down}}}}=cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}^{\mathrm{down}}}}\right)}^2}-{\displaystyle \frac{H_{\mathrm{wall}}^{\mathrm{down}}}{R_{\mathrm{wall}}^{\mathrm{down}}}}sin\left(i\right),$$

(3a)

$${\displaystyle \frac{\mathrm{i}{Y}_{\mathrm{down}}}{R_{\mathrm{wall}}^{\mathrm{down}}}}=-cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}^{\mathrm{down}}}}\right)}^2}-{\displaystyle \frac{H_{\mathrm{wall}}^{\mathrm{down}}}{R_{\mathrm{wall}}^{\mathrm{down}}}}sin\left(i\right).$$

(3b)

Ellipses ${\epsilon }_{\mathrm{up}}$ and ${\epsilon }_{\mathrm{down}}$ intersect at critical angles ${\theta }_{\mathrm{c}}$ y $\pi -{\theta }_{\mathrm{c}}$, where ${\theta }_{\mathrm{c}}$ is given by

$$sin{\theta }_{\mathrm{c}}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{H_{\mathrm{wall}}^{\mathrm{up}}}{R_{\mathrm{wall}}^{\mathrm{up}}}}+{\displaystyle \frac{H_{\mathrm{wall}}^{\mathrm{down}}}{R_{\mathrm{wall}}^{\mathrm{down}}}}\right)tan\left(i\right)=\delta .$$

(4)

Depending on the wall inclination angle i, there exist two possibilities to know whether both ellipses can intersect: if $\delta <1$ or not if $\delta >1$, as seen in Figure 3.

For the case $\delta <1$, the region $\mathcal{R}$ is composed by two sub-regions ${\mathcal{R}}_1$ and ${\mathcal{R}}_2$:

where $0<\theta <arcsin\left(\delta \right)$:

$${\mathcal{R}}_1=\left\{(X,Y):-X_0\le X\le X_0\wedge \mathrm{s}{Y}_{\mathrm{down}}\le Y\le \mathrm{s}{Y}_{\mathrm{up}}\right\},$$

(5a)

and where $arcsin\left(\delta \right)<\theta <\frac{\pi }{2}$:

$${\mathcal{R}}_2=\left\{(X,Y):\pm X_0\le X\le \pm R_{\mathrm{wall}}^{\mathrm{up}}\wedge \mathrm{i}{Y}_{\mathrm{up}}\le Y\le \mathrm{s}{Y}_{\mathrm{up}}\right\},$$

(5b)

that means $\mathcal{R}={\mathcal{R}}_1\cup {\mathcal{R}}_2$ (see Appendixes Appendix A.3 and Appendix A.4).

For the case $\delta >1$, the region $\mathcal{R}$ is defined as follows

$$\mathcal{R}=\left\{(X,Y):-R_{\mathrm{wall}}\le X\le R_{\mathrm{wall}}\wedge \mathrm{i}{Y}_{\mathrm{up}}\le Y\le \mathrm{s}{Y}_{\mathrm{up}}\right\}.$$

(6)

2.2. Vertical Walls

If in Definition 1 I set $R_{\mathrm{wall}}^{\mathrm{up}}=R_{\mathrm{wall}}^{\mathrm{down}}=R_{\mathrm{wall}}$, I obtain a two–dimensional vertical wall, as seen in Figure 4a. By rotating this line segment around z–axis, I generate a cylindrical ring which lies in the Euclidean space $(x,y,z)$, as seen in Figure 4b. This ring is a tri–dimensional cylindrical wall.

Following the same mathematical procedure as in the case of an inclined wall, I obtain that

$$\begin{array}{cc}\hfill {\epsilon }_{\mathrm{up}}:Y_{\mathrm{up}}& =\mathrm{s}{Y}_{\mathrm{up}}\cup \mathrm{i}{Y}_{\mathrm{up}},\hfill \end{array}$$

(7a)

$$\begin{array}{cc}\hfill {\epsilon }_{\mathrm{down}}:Y_{\mathrm{down}}& =\mathrm{s}{Y}_{\mathrm{down}}\cup \mathrm{i}{Y}_{\mathrm{down}},\hfill \end{array}$$

(7b)

where

$${\displaystyle \frac{\mathrm{s}{Y}_{\mathrm{up}}}{R_{\mathrm{wall}}}}=cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}}}\right)}^2}+{\displaystyle \frac{H_{\mathrm{wall}}}{R_{\mathrm{wall}}}}sin\left(i\right),$$

(8a)

$${\displaystyle \frac{\mathrm{i}{Y}_{\mathrm{up}}}{R_{\mathrm{wall}}}}=-cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}}}\right)}^2}+{\displaystyle \frac{H_{\mathrm{wall}}}{R_{\mathrm{wall}}}}sin\left(i\right),$$

(8b)

$${\displaystyle \frac{\mathrm{s}{Y}_{\mathrm{down}}}{R_{\mathrm{wall}}}}=cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}}}\right)}^2}-{\displaystyle \frac{H_{\mathrm{wall}}}{R_{\mathrm{wall}}}}sin\left(i\right),$$

(8c)

$${\displaystyle \frac{\mathrm{s}{Y}_{\mathrm{down}}}{R_{\mathrm{wall}}}}=cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}}}\right)}^2}-{\displaystyle \frac{H_{\mathrm{wall}}}{R_{\mathrm{wall}}}}sin\left(i\right).$$

(8d)

Both ellipses intersect at critical angles ${\theta }_{\mathrm{c}}$ y $\pi -{\theta }_{\mathrm{c}}$, where ${\theta }_{\mathrm{c}}$ is given by

$$sin{\theta }_{\mathrm{c}}={\displaystyle \frac{H_{\mathrm{wall}}}{R_{\mathrm{wall}}}}tan\left(i\right)=\delta .$$

(9)

3. The RHADaMAnTe Code

To create synthetic SEDs as arising from the inner curved wall of a gap or hole open by a planet in a protoplanetary disk, I have developed a computational code, written in the fortran 90 language, called rhadamante. This code is coupled to the artemise code because the geometry of the wall is required.

As I am interested in estimating the radiation reemitted by a tri–dimensional wall $\mathcal{W}$ projected on the plane of the sky, in this code, I firstly calculate the angle between the radial ray and the normal to the two–dimensional wall ${\mathcal{W}}_{{\pi }_0}$ for each incident radial radiation ray coming from the central star, as seen in Figure 5, by applying an algorithm also called RHADaMAnTe.

Then, I construct the tri–dimensional wall as the finite union of tri–dimensional conic rings obtained by rotating inclined line segments about the z–axis at different heights. (see Figure 6a and Figure 7).

Next, I calculate the surface projection on the plane of the sky of these rings, and then I calculate the radiation emitted by each of them by implementing some ideas from an algorithm developed for vertical walls [6]. Finally, I sum the contribution of the emission of all the projected rings to create a synthetic SED.

3.1. Geometry of the Radiation Reemitted by the Wall

rhadamante, acronym for Radial Geometry Algorithm for Calculating the Radiation Emitted by a Wall), it is a geometrical algorithm which at first calculates the angle between the stellar radiation along a radial ray ${\mathcal{R}}_k$ and the normal to the two–dimensional wall ${\mathcal{W}}_{{\pi }_0}$, as seen in Figure 5. Secondly, this algorithm discretizes the two–dimensional wall, which is not continuous, as seen in Figure 6a.

Let ${\theta }_{\mathit{k}}$ be the angle between the normal ${\eta }_{\mathit{k}}$ on the point $W_{\mathit{k}}$ belonging to the wall ${\mathcal{W}}_{{\pi }_0}$, and the stellar radiation ray ${\mathcal{R}}_{\mathit{k}}$, such that

$${\theta }_{\mathit{k}}={\theta }_{{\eta }_{\mathit{k}}}+{\phi }_k,$$

(10)

where ${\theta }_{{\eta }_{\mathit{k}}}$ is the minimal angle between the normal ${\eta }_{\mathit{k}}$ and the mid–plane (r–axis), and ${\phi }_k$ is the angle between the ray ${\mathcal{R}}_k$ and the mid–plane, as seen in Figure 5.

The angle ${\theta }_k$ is required to calculate the reemitted stellar radiation by the wall. Because of wall’s curvature and the radial geometry of the stellar radiation, each parcel of the wall does not absorb the total radiation, as it is in the case of vertical walls. In this case, each parcel absorbs only a fraction of the radiation which depends on the $cos\left({\theta }_k\right)$.

Let ${\ell }_{\mathit{k}}$ be the tangent line to the wall ${\mathcal{W}}_{{\pi }_0}$ at the point $W_{\mathit{k}}$ with a positive slope $m_{{\ell }_k}$. It follows that the inclination angle of such a line, measured from the r–axis, is ${\theta }_{{\ell }_k}={tan}^{-1}\left(m_{{\ell }_k}\right)$, where

$$m_{{\ell }_k}={\left.{\displaystyle \frac{\mathit{d}\mathcal{W}}{\mathit{d}\mathit{r}}}\right|}_{W_k}.$$

Physically, the wall ${\mathcal{W}}_{{\pi }_0}$ should be characterized by a mathematical continuous function. However, in this case, because of the numerical simulation, the wall is transfered, via the discretization process described in Section 3.2, into a discrete counterpart. So, as the points $W_k$ with $k=1,2,\cdots ,k_{max}$ defining ${\mathcal{W}}_{{\pi }_0}$ are close enough, it is possible to find an approximation of its derivative.

Consider the points $W_k$ and $W_{k+1}$ in the wall ${\mathcal{W}}_{{\pi }_0}$ to be connected along the segment line ${\mathcal{L}}_k$ (as seen in Figure 5), then the slope $m_{{\mathcal{L}}_k}$ of this line approximates to the derivative with respect to r of ${\mathcal{W}}_{{\pi }_0}$ at the point $W_k$, that is

$${\left.{\displaystyle \frac{\mathit{d}{\mathcal{W}}_{{\pi }_0}}{\mathit{d}\mathit{r}}}\right|}_{W_k}\approx m_{{\mathcal{L}}_k}.$$

Hence ${\theta }_{{\ell }_k}\approx {\theta }_{{\mathcal{L}}_k}={tan}^{-1}\left(m_{{\mathcal{L}}_k}\right)$.

Next, since the line ${\mathcal{L}}_k$ is almost perpendicular to the normal ${\eta }_{\mathit{k}}$, it follows

$${\theta }_{{\eta }_{\mathit{k}}}\approx \frac{\pi }{2}-{\theta }_{{\mathcal{L}}_k}.$$

(11)

Finally, since the star is located at the origin of coordinate system, it is easily to calculate the angle between the ray ${\mathcal{R}}_k$ and the mid–plane

$${\phi }_k={tan}^{-1}\left({\displaystyle \frac{W_k^{\left(z\right)}}{W_k^{\left(r\right)}}}\right),$$

(12)

where $W_k^{\left(r\right)}$ and $W_k^{\left(z\right)}$ are the r and z coordinates of the point $W_k$.

3.2. Discretization of the Two–Dimensional Wall

By applying the ARTeMiSE algorithm, I obtain a set of points $W_k=(r_k,z_k)$ with $k=1,2,\dots ,k_{max}$, defining the two–dimensional wall ${\mathcal{W}}_{{\pi }_0}$. It means that the wall is not a continuous curve. Then I discretize the wall as the finite union of infinitesimal inclined walls: I connect each couple of points $W_{k-1}$ and $W_k$ by inclined line segments ${\mathcal{W}}_{W_k}$ with height $2H_{\mathrm{wall}}^k$, as seen in Figure 6.

Consider an inclined line segment ${\mathcal{W}}_{W_k}$ with boundaries $W_{k-1}$ and $W_k$ in the wall ${\mathcal{W}}_{{\pi }_0}$, and let $Pm(W_{k-1},W_k)$ be the mid–point of ${\mathcal{W}}_{W_k}$. As I require that the vertical height of this line segment to be $2H_{\mathrm{wall}}^k$, I define

$$H_{\mathrm{wall}}^k:=W_k^{\left(z\right)}-P{m}^{\left(z\right)}(W_{k-1},W_k)=\frac{1}{2}(z_k-z_{k-1})\mathrm{if}k=2,3,\dots ,k_{max},$$

where $W_k^{\left(z\right)}$ and $P{m}^{\left(z\right)}(W_k,W_{k-1})$ are the z coordinates of $W_k$ and the mid–point between the points $W_k$ and $W_{k-1}$. See Figure 6.

If $k=1$, I construct a vertical line segment ${\mathcal{W}}_{W_1}$ with boundaries $W_0=(r_1,0)$ and $W_1$, and height $z_1$, which conects to the mid–plane.

3.3. Curved Wall

Let ${\mathcal{W}}_{{\pi }_0}={\bigcup }_k^{k_{max}}{\mathcal{W}}_{W_k}$ be a two–dimensional wall discretized by infinitesimal inclined line segments ${\mathcal{W}}_{W_k}$ for $k=2,\dots ,k_{max}$, and a vertical line segment ${\mathcal{W}}_{W_1}$.

By rotating each inclined line segment around z–axis, one generates a conic ring ${\mathcal{W}}_k$ with minimum radius $R_{k-1}$, maximum radius $R_k$ and total height $2H_{\mathrm{wall}}^k$. Whereas, by rotating the vertical line segment, a cylindrical ring with radius $R_1$ and total height $2H_{\mathrm{wall}}^1$, is obtained. It follows the tri–dimensional curved wall can be defined as the finite union of a cylindrical ring and several conic rings: $\mathcal{W}={\bigcup }_k^{k_{max}}{\mathcal{W}}_k$. See Figure 7.

3.3.1. Projection on the Plane of the Sky

The wall $\mathcal{W}$ has to be projected on the plane of the sky $(X,Y)$ to calculate the amount of visible surface. Therefore I consider the coordinate system $(X,Y,Z)$, where Z is along the line of sight, such that

$$X=x,$$

(13a)

$$Y=ycos\left(i\right)-zsin\left(i\right),$$

(13b)

$$Z=zcos\left(i\right)-ysin\left(i\right).$$

(13c)

Since I want to apply the same algorithms of projection described in Section 2.1 and Section 2.2, it is required to do a geometric translation of each ring ${\mathcal{W}}_k$ to a secondary coordinate system $(x^{\prime },y^{\prime },z^{\prime })$ such that the translated ring ${\mathcal{W}}_k^{\prime }$ is centered at the origin. Easily I can say that there exists a translation transformation ${\Theta }_1:(x,y,z)\to (x^{\prime },y^{\prime },z^{\prime })$:

$$x^{\prime }=x$$

(14a)

$$y^{\prime }=y$$

(14b)

$$z^{\prime }=z-\delta z_k,$$

(14c)

where $\delta z_k=z_k+\frac{1}{2}\Delta z_k$ is the displacement of the ring along z–axis, to be centered at the origin of the system $(x^{\prime },y^{\prime },z^{\prime })$, see Figure 8, with $z_k$ is the z coordinate of the point $W_k$ and $\Delta z_k=2H_{\mathrm{wall}}^k$.

Applying this translation, it is possible to use the coordinate system $(X^{\prime },Y^{\prime },Z^{\prime })$ to project the ring ${\mathcal{W}}_k^{\prime }$ on the plane $(X^{\prime },Y^{\prime })$:

$$X^{\prime }=x,$$

(15a)

$$Y^{\prime }=y^{\prime }cos\left(i\right)+z^{\prime }sin\left(i\right),$$

(15b)

$$Z^{\prime }=z^{\prime }cos\left(i\right)-y^{\prime }sin\left(i\right).$$

(15c)

From Equations (13), (14), and (15) it follows that there exists a translation transformation ${\Theta }_2:(X^{\prime },Y^{\prime },Z^{\prime })\to (X,Y,Z)$, such that:

$$X=X^{\prime },$$

(16a)

$$Y=Y^{\prime }+\delta z_{k}sin\left(i\right)$$

(16b)

$$Z=Z^{\prime }+\delta z_{k}cos\left(i\right),$$

(16c)

where $\delta z_{k}sin\left(i\right)$ is the projection on the Y–axis of the displacement of the ring along the z–axis, as seen in Figure 8.

Combining Equation (16) with Equations (2) and (3), it follows that for the k–th conic ring

$${\displaystyle \frac{\mathrm{s}{Y}_{\mathrm{up}}^{\left(k\right)}}{R_{\mathrm{wall},k}}}=cos\left(i\right)\sqrt{1-{\displaystyle \frac{X^2}{R_{\mathrm{wall},k}^2}}}+{\displaystyle \frac{H_{\mathrm{wall},k}+\delta z_k}{R_{\mathrm{wall},k}}}sin\left(i\right),$$

(17a)

$${\displaystyle \frac{\mathrm{i}{Y}_{\mathrm{up}}^{\left(k\right)}}{R_{\mathrm{wall},k}}}=-cos\left(i\right)\sqrt{1-{\displaystyle \frac{X^2}{R_{\mathrm{wall},k}^2}}}+{\displaystyle \frac{H_{\mathrm{wall},k}+\delta z_k}{R_{\mathrm{wall},k}}}sin\left(i\right).$$

(17b)

$${\displaystyle \frac{\mathrm{s}{Y}_{\mathrm{down}}^{\left(k\right)}}{R_{\mathrm{wall},k-1}}}=cos\left(i\right)\sqrt{1-{\displaystyle \frac{X^2}{R_{\mathrm{wall},k-1}^2}}}-{\displaystyle \frac{H_{\mathrm{wall},k-1}-\delta z_{k-1}}{R_{\mathrm{wall},k-1}}}sin\left(i\right),$$

(17c)

$${\displaystyle \frac{\mathrm{i}{Y}_{\mathrm{down}}^{\left(k\right)}}{R_{\mathrm{wall},k-1}}}=-cos\left(i\right)\sqrt{1-{\displaystyle \frac{X^2}{R_{\mathrm{wall},k-1}^2}}}-{\displaystyle \frac{H_{\mathrm{wall},k-1}-\delta z_{k-1}}{R_{\mathrm{wall},k-1}}}sin\left(i\right).$$

(17d)

Whereas, for the cylindrical ring, by combining Equation (16) with Equation (8), it follows

$${\displaystyle \frac{\mathrm{s}{Y}_{\mathrm{up}}^{\left(1\right)}}{R_{\mathrm{wall},k}}}=cos\left(i\right)\sqrt{1-{\displaystyle \frac{X^2}{R_{\mathrm{wall},k}^2}}}+{\displaystyle \frac{H_{\mathrm{wall},k}+\delta z_k}{R_{\mathrm{wall},k}}}sin\left(i\right),$$

(18a)

$${\displaystyle \frac{\mathrm{i}{Y}_{\mathrm{up}}^{\left(1\right)}}{R_{\mathrm{wall},k}}}=-cos\left(i\right)\sqrt{1-{\displaystyle \frac{X^2}{R_{\mathrm{wall},k}^2}}}+{\displaystyle \frac{H_{\mathrm{wall},k}+\delta z_k}{R_{\mathrm{wall},k}}}sin\left(i\right).$$

(18b)

$${\displaystyle \frac{\mathrm{s}{Y}_{\mathrm{down}}^{\left(1\right)}}{R_{\mathrm{wall},k}}}=cos\left(i\right)\sqrt{1-{\displaystyle \frac{X^2}{R_{\mathrm{wall},k}^2}}}-{\displaystyle \frac{H_{\mathrm{wall},k}-\delta z_k}{R_{\mathrm{wall},k}}}sin\left(i\right),$$

(18c)

$${\displaystyle \frac{\mathrm{i}{Y}_{\mathrm{down}}^{\left(1\right)}}{R_{\mathrm{wall},k}}}=-cos\left(i\right)\sqrt{1-{\displaystyle \frac{X^2}{R_{\mathrm{wall},k}^2}}}-{\displaystyle \frac{H_{\mathrm{wall},k}-\delta z_k}{R_{\mathrm{wall},k}}}sin\left(i\right).$$

(18d)

3.4. Emission of the Wall

To calculate the emission or emergent flux $F_{\nu }$ of the visible wall projected on the plane of the sky, I multiply the total emergent intensity $I_{\nu }$ by the solid angle ${\Omega }_{\mathrm{wall}}$ of the visible surface of the wall, whose geometry has been described in detail in Section 3.3.

For each element in the visible surface of the projected wall, the thermal emergent intensity, approximated as isotropic, is given by

$$I_{\nu }\approx {\int }_0^{\infty }{B}_{\nu }\left[T_{\mathrm{d}}\left({\tau }_{\mathrm{d}}\right)\right]exp\left(-{\tau }_{\nu }\right)d{\tau }_{\nu }$$

(19)

[see [6] for derivation], where $B_{\nu }$ is the Planck function, ${\tau }_{\mathrm{d}}$ is the total mean optical depth at the disk frequency band, and ${\tau }_{\nu }={\tau }_{\mathrm{d}}({\kappa }_{\nu }/{\chi }_{\mathrm{d}})$, with opacity ${\kappa }_{\nu }$.

The wall temperature $T_{\mathrm{d}}$ is a function of the optical depth of the disk, and it is calculated as follows [see [6] for derivation]

$$T_{\mathrm{d}}^4\left({\tau }_{\mathrm{d}}\right)=\alpha {\displaystyle \frac{F_0}{4{\sigma }_{\mathrm{R}}}}\left[C_1^{\prime }+C_2^{\prime }exp(-q{\tau }_{\mathrm{d}})+C_3^{\prime }exp(-\beta q{\tau }_{\mathrm{d}})\right],$$

(20)

where

$$C_1^{\prime }=(1+C_1)\left(2+{\displaystyle \frac{3}{q}}\right)+C_2\left(2+{\displaystyle \frac{3}{\beta q}}\right),$$

(21a)

$$C_2^{\prime }=(1+C_1)\left({\displaystyle \frac{q{\chi }_{\mathrm{d}}}{{\kappa }_{\mathrm{d}}}}-{\displaystyle \frac{3}{q}}\right),$$

(21b)

$$C_3^{\prime }=C_2\beta \left({\displaystyle \frac{q{\chi }_{\mathrm{d}}}{{\kappa }_{\mathrm{d}}}}-{\displaystyle \frac{3}{q{\beta }^2}}\right),$$

(21c)

and

$$C_1=-{\displaystyle \frac{3w}{1-{\beta }^2}},$$

(22a)

$$C_2={\displaystyle \frac{5w}{\beta \left[(1+\frac{2\beta }{1}\right](1-{\beta }^2)}},$$

(22b)

with $\alpha =1-w$, $\beta =\sqrt{3\alpha }$, and $w={\sigma }_{\mathrm{s}}/{\chi }_{\mathrm{s}}$ is the mean albedo to the stellar radiation and $F_0=L_{★}/4\pi R_{\mathrm{wall}}^2$, where $L_{★}$ is the stellar luminosity.

At a distance d from the observer, the total solid angle is given by

$${\Omega }_{\mathrm{wall}}=\{\begin{array}{c}cos\left(i\right){\left({\displaystyle \frac{R_{\mathrm{wall}}^{\mathrm{up}}+R_{\mathrm{wall}}^{\mathrm{down}}}{d}}\right)}^2\left[\delta \sqrt{1-{\delta }^2}+arcsin\left(\delta \right)\right],\mathrm{if}\delta <1,\hfill \\ \pi cos\left(i\right){\left({\displaystyle \frac{R_{\mathrm{wall}}^{\mathrm{up}}}{d}}\right)}^2,\mathrm{if}\delta \ge 1,\hfill \end{array}$$

(23)

with

$$\delta ={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{H_{\mathrm{wall}}^{\mathrm{up}}}{R_{\mathrm{wall}}^{\mathrm{up}}}}+{\displaystyle \frac{H_{\mathrm{wall}}^{\mathrm{down}}}{R_{\mathrm{wall}}^{\mathrm{down}}}}\right)tan\left(i\right),$$

for conic rings, and

$${\Omega }_{\mathrm{wall}}=\{\begin{array}{c}2cos\left(i\right){\left({\displaystyle \frac{R_{\mathrm{wall}}}{d}}\right)}^2\left[\delta \sqrt{1-{\delta }^2}+arcsin\left(\delta \right)\right],\mathrm{if}\delta <1,\hfill \\ \pi cos\left(i\right){\left({\displaystyle \frac{R_{\mathrm{wall}}}{d}}\right)}^2,\mathrm{if}\delta \ge 1,\hfill \end{array}$$

(24)

with

$$\delta ={\displaystyle \frac{H_{\mathrm{wall}}}{R_{\mathrm{wall}}}}tan\left(i\right),$$

for cylindrical rings.

3.4.1. Rosseland Mean Opacity

Equation (19) requires the calculation of the opacity ${\kappa }_{\nu }$. This dominant opacity depends on the chemical composition, pressure and temperature of the gas, as well as the frequency $\nu$ of the incident light. This is a complex endeavour. The problem can be simplified by using a mean opacity averaged over all frequencies, so that only the dependence on the gas physical properties remains. In the current work, I use the Rosseland mean opacity, defined as

$${\displaystyle \frac{1}{\langle {\kappa }^R\rangle }}:={\displaystyle \frac{{\displaystyle {\int }_0^{\infty }\frac{1}{{\kappa }_{\nu }}\frac{\partial B_{\nu }}{\partial T}d\nu }}{{\displaystyle {\int }_0^{\infty }\frac{\partial B_{\nu }}{\partial T}d\nu }}},$$

(25)

where $B_{\nu }\left(T\right)$ is the Planck’s function, and T is the disk temperature [32].

To calculate the total Rosseland mean opacity ${\kappa }^{\mathrm{R}}$, I consider that all the dust grains species exist and the mixture of dust grains is made of small and big grains. Using the previous assumptions I calculate the total Rosseland mean opacity as follows:

$${\kappa }_{\mathrm{R}}(x,z)={\zeta }_{\mathrm{small}}(x,z){\kappa }_{{\mathrm{R}}_{\mathrm{small}}}+{\zeta }_{\mathrm{big}}(x,z){\kappa }_{{\mathrm{R}}_{\mathrm{big}}},$$

(26)

where ${\kappa }_{{\mathrm{R}}_{\mathrm{small}}}$ and ${\kappa }_{{\mathrm{R}}_{\mathrm{big}}}$ are the Rosseland mean opacities associated to the small dust grain size distribution and big dust grain size distribution, respectively. And ${\zeta }_{\mathrm{small}}$ and ${\zeta }_{\mathrm{big}}$ represent the abundances (dust-to-gas mass ratio) of the small and big grains, respectively:

$${\zeta }_{\mathrm{small}}(x,z)=\frac{1}{2}{\zeta }_{\mathrm{small},0}\left\{1-tanh\left[k\left(1-{\displaystyle \frac{z}{\delta H}}\right)\right]\right\},$$

(27a)

$${\zeta }_{\mathrm{big}}(x,z)=\frac{1}{2}{\zeta }_{\mathrm{big},0}\left\{1+tanh\left[k\left(1-{\displaystyle \frac{z}{\delta H}}\right)\right]\right\},$$

(27b)

here $\delta H$ represents a small fraction of the scale height of the disk, and k is a factor which defines a smooth transition between small and big grains population [33].

The monochromatic opacity ${\kappa }_{\nu }$ in Equation (25) depends on the dust species in the mixture and their physical and chemical properties, such that it is calculated as the sum of the monochromatic opacity of each grain species:

$${\kappa }_{\nu }=\sum _q^{q_{max}}{\kappa }_{\nu }^q(a_{min}^q,a_{max}^q,{\sigma }^q,{\eta }^q),$$

(28)

where $a_{min}^q$ and $a_{max}^q$ are the sizes of the small and big grains, and ${\sigma }^q$ and ${\eta }^q$ are the abundance and refraction index of the species. Here q is running over the name of the species (e.g. silicates, organics, amorphous carbon, ice and troilite) in the dust composition of the disk. I calculate the monochromatic opacities using the Mie theory by implementing some modified routines of a code developed in [6].

Summarizing, to calculate the emergent flux emitted by the curved wall, I have developed a computational code called rhadamante. This code is based on the geometry of the wall calculated by the artemise code. In Figure 9 I show a flowchart of our code. For some tests, see Appendix B.

4. Results: An Implementation to the Stellar System LkCa 15

In this section, I present a model of the truncated dusty disk of the T Tauri star LkCa 15 that accounts for the Spitzer Infrared Spectrograph observations. I have modeled the mid–infrared spectral energy distribution from 5 and 40$\mu \mathrm{m}$) as arising from the inner curved wall of the outer disk. In this model a $10M_{♃}$ mass planet is the responsible of the wall curvature. The free dust hole has a radius of ∼$53\mathrm{AU}$ along the mid-plane. The wall has a half-height of ∼$12\mathrm{AU}$ and it is illuminated at normal incidence by the central star, but it also is shadowed because of the presence of an internal optically thick disk.

4.1. Simulation: Planet–Disk Interaction

As I am interested in characterizing the geometry of the wall of the disk gap, in the LkCa 15 system, I need to analyze the vertical structure of the disk. Assuming the gap was opened by an embedded planet, I use the fargo–3d hydro–dynamical code to launch two numerical simulations of the disk–planet interaction until the 500th. The only difference among these simulations is the size in resolution ($N_{\mathrm{X}}\times N_{\mathrm{Y}}\times N_{\mathrm{Z}}$). The low resolution of $50\times 50\times 30$ was used to find quickly the orbit where the system reaches a quasi–stationary state. The medium resolution of $250\times 250\times 100$ was used to get a better approximation of the wall. According to fargo–3d requirements, in Table 1 I show some parameters for simulations.

In Figure 10, I present the tri–dimensional structure of the 100th orbit of the LkCa 15 disk simulation. It is here where the system reached a stationary-state.

In Figure 11, it can be observed an optically thick inner disk and the gap completly opened. The planet is located at the point $(x,y)=(1,0)$ as seen in Figure 11.

4.2. Dust Grain

The optical depth of the disk depends on the disk material opacity, I assume that the disk is a mixture of grains composed of silicates with mass fraction ${\zeta }_{\mathrm{sil}}=0.0034$, organics with ${\zeta }_{\mathrm{org}}=0.0041$, and troilite with ${\zeta }_{\mathrm{troi}}=8.0\times {10}^{-4}$, consistent with the model proposed in [34]. The grains are assumed to be spheres, which obey the standard MRN grain size distribution $n\left(a\right)\sim$$a^{3.5}$ [7].

I consider two grain populations: small grains between minimum radius $a_{min}=0.005\mu \mathrm{m}$ and maximum radius $a_{max}=0.25\mu \mathrm{m}$, and big grains between minimum radius $a_{min}=0.005\mu \mathrm{m}$ and maximum radius $a_{max}=1000\mu \mathrm{m}$. I consider a smooth transition between both dust populations, see Equation (27) in Section 3.4.1, where I set $k=20$, $\delta =0.1$ and ${\zeta }_{\mathrm{small},0}=0.5$ and ${\zeta }_{\mathrm{big},0}=6.8$. I use optical constants for silicates from [35,36] and [37], for the organics from [34] and for troilite from [34] and [38]. I also take into account the sublimation temperature of the grain species in the mixture $T_{\mathrm{sil}}=1400\mathrm{K}$, $T_{\mathrm{org}}=425\mathrm{K}$, and $T_{\mathrm{troi}}=680\mathrm{K}$.

For the composition of the silicate dust grains, I considered several possibilities as showed in [39]: amorphous Mg–Fe glassy olivine (${\mathrm{Mg}}_{2x}{\mathrm{Fe}}_{2-2x}{\mathrm{SiO}}_4$) and glassy pyroxene (${\mathrm{Mg}}_x{\mathrm{Fe}}_{1-x}{\mathrm{SiO}}_3$), where the subscripts to each element represents the fraction of the silicate composed of that element, with optical constants from [36].

4.3. The Vertical Geometry of the Wall

I found that the $10M_{♃}$ mass planet candidate, when is located at $32.3\mathrm{AU}$ from the central star, opens a gap around the young transitional disk host LkCa 15. The artemise code was implemented to analize the simulation data. The radii of the wall along the mid–plane of the disk and the heights of the wall have a deppendence on the chemical composition of the silicate grains as showed in Table 2. In Figure 12, I show the geometry of a wall where the dust garin disk composition has glassy olivine (silicate) with 50%Fe and 50% Mg.

The location of the planet candidate is not consistent with the observations [e.g. [25], which suggest that the possible massive planets LkCa 15 b and LkCa 15c are located at $14.7\pm 2.1\mathrm{AU}$ and $18.6\pm 2.5\mathrm{AU}$, respectively, along the semimajor axis. However, the radii of the wall along the mid–plane $R_{\mathrm{wall}}^{\mathrm{mid}-\mathrm{plane}}$ are similar to those measured in [40,41], ∼$50\mathrm{AU}$, and [22],∼$56\mathrm{AU}$.

4.4. SED of the Wall

I model LkCa 15 as a central star with the properties described in Section 1.2, surrounded by an optically thick inner disk (as showed in Figure 11) and and outer disk truncated at ∼$120\mathrm{AU}$. I consider the gap has a curved wall at differentent locations and heights according to Table 2. In the models, I consider that LkCa 15 is at $140\mathrm{pc}$ from Earth in the Taurus-Auriga star forming region [21] and the disk inclination is $i={50}^{\circ }$ [25]. A representation of the model is showed in Figure 13.

As the inner disk casts an umbra over the wall of the outer disk, as seen in Figure 13, I have to remove it from the SED of the outer disk wall. In order to find $z_{\mathrm{umb}}$, I implement some improved routines developed in [42] for curved sublimation walls. This code uses opacities to calculate the shape of the wall and assumes that the stellar rays are parallel to the mid–plane. I found that the wall of the inner disk starts at ∼$0.04932\mathrm{AU}$, from the central star, and runs until ∼$0.18407\mathrm{AU}$ where it reaches ∼$2.3550\times {10}^{-2}\mathrm{AU}$ in height. And the temperature of the sublimation wall decreases with radius and it ranges from $1014.3\mathrm{K}$ to $1443.7\mathrm{K}$.

In a first approximation, assuming the star as a point, the sublimation wall produces only an umbra over the wall of the outer disk (see Figure 13). To calculate the size $z_{\mathrm{umb}}$ of this umbra, theangle ${\theta }_{\mathrm{shw}}$ subtended by the height $z_{\mathrm{wall}}^{\mathrm{inner}}$ of the sublimation wall is needed. In addition, for some points $(r,z)$ in the outer wall, I calculate the angle ${\theta }_{\mathrm{shw}}^{\prime }=arctan(z/r)$ until it reaches the value ${\theta }_{\mathrm{shw}}$.

Applying the previous algorithm to the geometry of the wall (see Figure 12), I found that the umbra produced by the sublimation wall of the inner disk onto the wall of the outer disk is $7.44\mathrm{AU}$ in height (above and below the disk mid–plane). It means the contribution to the SED of the outer wall comes from a region of the wall from ∼$58.84\mathrm{AU}$ to ∼$68.7\mathrm{AU}$ along the radial direction, and from $7.44\mathrm{AU}$ to $12.22\mathrm{AU}$ along the vertical direction. In Figure 14a I show the surface area of the whole curved wall, and in Figure 14b I show the surface area considered the umbra cast by the inner disk.

Synthetic SEDs of the wall of the outer disk, where the dust consists of grains of glassy pyroxene with different concentration of Fe and Mg are showed in Figure 15. I performed a chi–square test for each model, to examine whether the synthetic SED fits to the Spitzer IRS SED. I found that no model is either capable to fit the observed SED (${\chi }^2$∼300) nor reproduce the silicate peak at ∼$10\mu \mathrm{m}$. However, for glassy pyroxene with 60% Fe and 40% Mg (see Figure 15a) it seems the silicate peaks tries to appear at ∼$10\mu \mathrm{m}$, this lead us to think that a lesser concentration of Mg in the pyroxene composition would produce the silicate feature. Unfortunately, there is not available, in the literature, the optical constants needed to calculate the opacities, for such chemical concentrations.

In Figure 16, I show the synthetic SEDs of the wall of the outer disk, where the dust mixture consists of grains of glassy olivine with different concentration of Fe and Mg. I found that none of these configurations is capable to fit the observed SED (${\chi }^2>200$). However, in both cases a silicate peak appears at ∼$10\mu \mathrm{m}$. A concentration of 50% Fe and 50% Mg produces the best fit (see Figure 16b).

Previous results lead us to say that the SED of LkCa 15 is not dominated by the contribution of the curved wall of the outer disk in the mid–infrared. However, when olivine grains with a concentration of 50% Fe and 50% Mg or 80% Fe and 20% Mg are in the dust mixture, a silicate feature appears at ∼$10\mu \mathrm{m}$.

4.5. SED of the System

Considering the SED contribution of the inner sublimation wall and the star, in addition to the SED contribution of the wall of the outer disk, I present a more complete model of the stellar system LkCa 15. For the stellar SED, I have used a Kurucz atmosphere model with $log[Z/H]=0.5$ and $log\left[g\right]=4.0$. Table 3 lists the parameters for this model. In Figure 17, I show the contribution of the star, the inner sublimation wall and the wall of the outer disk to the total synthetic SED.

Only for some wavelengths in the mid–infrared, ∼$15.0\mu \mathrm{m}<\lambda <$∼$20.0\mu \mathrm{m}$, the sythetic SED fit to LkCa 15 observed SED, I estimated ${\chi }^2$∼$0.45$. For all the other wavelengths in the field of view of the Spitzer IRS (5.217–37.86 $\mu$m), the synthetic SED is below the observed SED. For wavelengths ∼$5.217\mu \mathrm{m}<\lambda \le$∼$8.0\mu \mathrm{m}$ the difference is not very high (${\chi }^2$∼$3.36$); however, for wavelengths ∼$8.0\mu \mathrm{m}<\lambda \le$∼$15.0\mu \mathrm{m}$ (with ${\chi }^2$∼$17.72$) and ∼$20.0\mu \mathrm{m}<\lambda \le$∼$37.86\mu \mathrm{m}$ (with ${\chi }^2$∼$19.08$) this difference becomes significant (see Figure 18).

I can suggests that the inner sublimation wall and the stellar photo-sphere cannot account for the significant near-infrared excess in LkCa 15. It means the SED model also requires of optically thin dust inside the gap to explain the excess and to produce the $10\mu \mathrm{m}$ silicate feature as showed in [43]. Similarly, as the wall of the outer disk cannot account for the excess in the mid and long–infrared, I asumme that this SED model also requires of the contribution of the outer disk.

5. Discussion: Vertical Wall SED vs. Curved Wall SED

Wall SED models of LkCa 15 (and many other pre-TDs and TDs) are based on vertical walls [e.g. [23,43]]. rhadamante code is able to construct wall SEDs based on this geometry. Here I present some wall SED models of LkCa 15 considering vertical walls to compare with the best fit curved wall SED model (as seen in Figure 19).

In these vertical wall SED models, it is also cosidered that the umbra cast by the inner disk is $6.73\mathrm{AU}$ (measured from the disk mid–plane to up). The height of the wall is $12\mathrm{AU}$ (measured from the disk mid–plane to up). The size and composition of the dust grains remain the same as described for the best fit curved wall SED model.

The vertical walls were located at 30, $58.11$, $60.85$ and $68.7\mathrm{AU}$. rhadamante estimated the area of the visible surface projected on the plane of the sky of these vertical walls, $A_{\mathrm{wall}}^{\mathrm{vert}}$, and the curved wall, $A_{\mathrm{wall}}^{\mathrm{cur}}$. I found that $A_{\mathrm{wall}}^{\mathrm{vert}}>A_{\mathrm{wall}}^{\mathrm{cur}}$, if the umbra of the inner disk is considered. While if the emission of the whole wall is considered, $A_{\mathrm{wall}}^{\mathrm{vert}}<A_{\mathrm{wall}}^{\mathrm{cur}}$. Wall temperatures also were estimated, I found $27.63$, $26.41$, $25.80$ and $24.29\mathrm{K}$, respectively. It means that the temperature of a vertical wall, $T_{\mathrm{wall}}$, decreases if its radius, $R_{\mathrm{wall}}$, increases. Furthermore, the radiation emitted by a vertical wall, $F_{\lambda }$, also decreases as $R_{\mathrm{wall}}$ increases, as seen in Figure 19.

Finaly, I compared the radiation emitted by vertical walls, $F_{\lambda }^{\mathrm{vert}}$, with the radiation emitted by the curved wall, $F_{\lambda }^{\mathrm{curved}}$ (see Figure 19). I found that $F_{\lambda }^{\mathrm{vert}}>F_{\lambda }^{\mathrm{curved}}$ for wavelengths between 5 and $1000\mu \mathrm{m}$. The difference between fluxes becomes significant, about one order of magnitude, for $9\mu \mathrm{m}<\lambda <35\mu \mathrm{m}$.

This infrared excess arises, in part, from the angle, $\theta$, between the radiation ray and the normal to the wall (see Section 3.1), because $F_{\lambda }\propto cos\theta$. For vertical walls, $cos\left(\theta \right)\sim 1$ for all the radiation rays hitting the wall, because the normal to the wall is always parallel to the disk mid–plane. Whereas for curved walls $0\le cos\left(\theta \right)\le 1$, because of the wall curvature. In addition, for some models, some regions of the curved wall are farther from the central star than the vertical one.

The above two facts induce a lower exposure of the curved wall to the host star radiation, which derives in much less radiative heating of the wall, and, consequently, in the significantly lower radiative infrared cooling flux. For this reason, the one order of magnitude in lower infrared emisison in the curved wall model is as significant as correct by very basic physical considerations: the height of the curved wall and the umbra cast by the inner disk onto the outer disk.

It is worth mentioning that the geometry, location and height of the curved wall arise from a physical mechanism considering the opacities and chemical composition of the disk, where the disk is the result of a three–dimensional hydro–dynamical simulation, whereas for the vertical walls there is not any physics to choice these parameters. Although the curved wall SED model is the best choice to compare with the Spitzer IRS SED of LkCa 15, as it is showed in Figure 19, to rise this synthetic flux to the right level, it is needed, in adittion, the contribution of the central star, the inner disk, the outer disk and a region of optically thin dust inside the gap, as discussed in Section 4.4.

6. Conclusions

The computational code rhadamante was developed in order to calculate synthetic SEDs of protoplanetry disks. As a initial parameter, the code requires the geometry of the wall coming from hydro-dynamical three-dimensional simulations of the planet-disk interaction and dust grain properties. This code is useful to explain the observed SED of young stellar systems in transition stage. It would lead to unveil the structure of the system, such as, the inner disk, the gap, and the outer disk, even the location and mass of the embedded planet responsible of the gap opening.

From the implementation of this code to the pre-transitional disk LkCa 15, it can be concluded that:

All models of the SED consisting only of a curved wall, using different concentrations of Fe and Mg for the silicate (pyroxene and olivine) grains, suggest that LkCa 15’s Spitzer IRS SED cannot be accounted by the emission of a curved wall. Chi-square tests indicate the models are not good (${\chi }^2$∼${10}^2$). In addition, for models with a dust mixture containing glassy amorphous olivine grains, ${\mathrm{Mg}}_{0.2}{\mathrm{Fe}}_{0.2}{\mathrm{SiO}}_4$ or ${\mathrm{Mg}}_{0.5}{\mathrm{Fe}}_{0.5}{\mathrm{SiO}}_4$, a $10\mu \mathrm{m}$ silicate feature can be observed. However, the intensity of this silicate emission to that measured in the Spitzer IRS SED.

Including to the wall SED model, the contribution of the stellar photo-sphere and the sublimation wall of the inner disk, considering a dust mixture with amorphous glassy olivine with 50% Fe and 50% Mg and a small amount of organic and troilite grains, the observed and synthetic SEDs fit better (${\chi }^2$∼45). However, this model cannot fit the Spitzer IRS SED. Only for a small band in the mid-infrared, ∼$15.5\mu \mathrm{m}<\lambda <$∼$18.0\mu \mathrm{m}$, the fit is good (${\chi }^2$∼$0.076$).

Several limitations exist in the current model: The contribution of the inner disk was not considered, neither the contribution of an optically thin region inside the gap [e.g. [23,43], nor the optically thick outer disk. The inner disk might contribute to the SED at wavelengths in the near infrared. The optically thin region might explain the silicate feature of the Spitzer IRS SED at $10\mu \mathrm{m}$. While the outer disk might contribute to the SED at wavelengths longer than $18\mu \mathrm{m}$.

Vertical wall SED models, via rhadamante code, show a difference of one order of magnitude in the flux, $F_{\lambda }^{\mathrm{vert}}$, compared to the curved wall SED model, $F_{\lambda }^{\mathrm{curved}}$, for wavelengths from 8 to $35\mu \mathrm{m}$. This difference arises from dependency of the flux, $F_{\lambda }$, on the cosine of the angle, $\theta$, between the stellar radiation ray and the normal to the wall. This lower exposure of the curved wall to the stellar radiation results in much less radiative heating of the wall, and, consequently, in the significantly lower radiative infrared cooling flux.

The synthetic SED of a curved wall, estimated by rhadamante code, includes physical and chemical mechanisms absent in the estimation of the SED of a vertical wall. That is, a curved wall SED is better to fit the Spitzer IRS SED of LkCa 15 or any other protoplanetaray disk.

Appendix A.1. Construction of a Vertical Wall

Consider an inclined segment line as showed in Figure 1a. Circles ${\mathcal{C}}_{\mathrm{up}}$ and ${\mathcal{C}}_{\mathrm{down}}$, with radii $R_{\mathrm{wall}}^{\mathrm{up}}$ and $R_{\mathrm{wall}}^{\mathrm{down}}$ and height $H_{\mathrm{wall}}^{\mathrm{up}}$ and $H_{\mathrm{wall}}^{\mathrm{down}}$, respectively, are generated by rotating such segment line around z–axis. These circles define the superior and inferior boundaries of a tri–dimensional wall $\mathcal{W}$. The projection on the $XY$–plane of these circles, given an inclination angle i, generates the ellipses ${\epsilon }_{\mathrm{up}}$ and ${\epsilon }_{\mathrm{down}}$, as seen in Figure A1a. Ellipse ${\epsilon }_{\mathrm{up}}$ is centered at $(h,k)$ with semi–major and semi–minor axes a and b, respectively. Whereas, ellipse ${\epsilon }_{\mathrm{down}}$ is centered at $(h^{\prime },k^{\prime })$ with semi–major and semi–minor axes $a^{\prime }$ and $b^{\prime }$, respectively.

Without lost of generality I focus on the construction of ellipse ${\epsilon }_{\mathrm{up}}$. As the coordinate system $XY$ is centered on the star, it follows $h=0$, and the Y coordinate of the ellipse is the projection of the height $H_{\mathrm{wall}}^{\mathrm{up}}$ of the circle ${\mathcal{C}}_{\mathrm{up}}$ on the Y-axis, as seen in Figure A1b, that is, $k=H_{\mathrm{wall}}^{\mathrm{up}}sin\left(i\right)$. The semi–major axis is $a=R_{\mathrm{wall}}^{\mathrm{up}}$, and the semi–minor axis is the projection of the radius of the circle ${\mathcal{C}}_{\mathrm{up}}$ on the Y-axis, as seen in Figure A1b, that is, $b=R_{\mathrm{wall}}^{\mathrm{up}}cos\left(i\right)$.

Figure A1. Geometry of ellipses ${\epsilon }_{\mathrm{up}}$ and ${\epsilon }_{\mathrm{down}}$ defining the projection of the tri–dimensional wall $\mathcal{W}$ on the $XY$–plane.

Figure A1. Geometry of ellipses ${\epsilon }_{\mathrm{up}}$ and ${\epsilon }_{\mathrm{down}}$ defining the projection of the tri–dimensional wall $\mathcal{W}$ on the $XY$–plane.

So the full form of the equation of ellipse ${\epsilon }_{\mathrm{up}}$ is:

$${\left[{\displaystyle \frac{X}{R_{\mathrm{wall}}^{\mathrm{up}}}}\right]}^2+{\left[{\displaystyle \frac{Y-H_{\mathrm{wall}}^{\mathrm{up}}sin\left(i\right)}{R_{\mathrm{wall}}^{\mathrm{up}}cos\left(i\right)}}\right]}^2=1.$$

(A1)

Similarly, for ellipse ${\epsilon }_{\mathrm{down}}$,

$${\left[{\displaystyle \frac{X}{R_{\mathrm{wall}}^{\mathrm{down}}}}\right]}^2+{\left[{\displaystyle \frac{Y+H_{\mathrm{wall}}^{\mathrm{down}}sin\left(i\right)}{R_{\mathrm{wall}}^{\mathrm{down}}cos\left(i\right)}}\right]}^2=1.$$

(A2)

Hence, I define

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{s}{Y}_{\mathrm{up}}}{R_{\mathrm{wall}}^{\mathrm{up}}}}& =cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}^{\mathrm{up}}}}\right)}^2}+{\displaystyle \frac{H_{\mathrm{wall}}^{\mathrm{up}}}{R_{\mathrm{wall}}^{\mathrm{up}}}}sin\left(i\right),\hfill \end{array}$$

(A3a)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{i}{Y}_{\mathrm{up}}}{R_{\mathrm{wall}}^{\mathrm{up}}}}& =-cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}^{\mathrm{up}}}}\right)}^2}+{\displaystyle \frac{H_{\mathrm{wall}}^{\mathrm{up}}}{R_{\mathrm{wall}}^{\mathrm{up}}}}sin\left(i\right),\hfill \end{array}$$

(A3b)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{s}{Y}_{\mathrm{down}}}{R_{\mathrm{wall}}^{\mathrm{down}}}}& =cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}^{\mathrm{down}}}}\right)}^2}-{\displaystyle \frac{H_{\mathrm{wall}}^{\mathrm{down}}}{R_{\mathrm{wall}}^{\mathrm{down}}}}sin\left(i\right),\hfill \end{array}$$

(A3c)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{i}{Y}_{\mathrm{down}}}{R_{\mathrm{wall}}^{\mathrm{down}}}}& =-cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}^{\mathrm{down}}}}\right)}^2}-{\displaystyle \frac{H_{\mathrm{wall}}^{\mathrm{down}}}{R_{\mathrm{wall}}^{\mathrm{down}}}}sin\left(i\right),\hfill \end{array}$$

(A3d)

such that

$$\begin{array}{cc}\hfill {\epsilon }_{\mathrm{up}}:Y_{\mathrm{up}}& =\mathrm{s}{Y}_{\mathrm{up}}\cup \mathrm{i}{Y}_{\mathrm{up}},\hfill \end{array}$$

(A4a)

$$\begin{array}{cc}\hfill {\epsilon }_{\mathrm{down}}:Y_{\mathrm{down}}& =\mathrm{s}{Y}_{\mathrm{down}}\cup \mathrm{i}{Y}_{\mathrm{down}}.\hfill \end{array}$$

(A4b)

Appendix A.2. Area Between Two Curves

Definition A1.

The area between the curves $f\left(x\right)$ and $g\left(x\right)$ and the ordinates $x=a$ and $x=b$ is given by

$$A={\int }_a^b\left[f\left(x\right)-g\left(x\right)\right]dx,$$

(A5)

if and only if $f\left(x\right)\ge g\left(x\right)$$\forall x\in [a,b]$.

Definition A2.

Let $f:\mathbb{R}\to [0,1]$ be a continuous function defined by $x\mapsto \sqrt{1-{(x/a)}^2}$. An anti-derivative or primitive function of f is

$$F\left(x\right)=\int \sqrt{1-{\left({\displaystyle \frac{x}{a}}\right)}^2}dx=\frac{1}{2}x\sqrt{1-{\left({\displaystyle \frac{x}{a}}\right)}^2}+\frac{1}{2}a{sin}^{-1}\left({\displaystyle \frac{x}{a}}\right)+C,$$

(A6)

such that $F^{\prime }\left(x\right)=f\left(x\right)$$\forall x\in [0,1]$.

Appendix A.3. Theoretical area of the projected vertical wall: The whole wall

The following analysis, on the area of a visible wall (projected on the plane of the sky), focuses on the case corresponding to $\delta <1$ (the star is visible), because of LkCa 15 inclination angle ($i={50}^{\circ }$).

The area of the projection of a whole vertical wall on the plane of the sky is

$$A=A_1+2A_2$$

(A7)

as seen in Figure (Figure A2). Where, according to Definition A1,

$$A_1={\int }_{-X_0}^{X_0}(\mathrm{s}{Y}_{\mathrm{up}}-\mathrm{s}{Y}_{\mathrm{down}})dX,$$

(A8a)

$$A_2={\int }_{X_0}^{R_{\mathrm{wall}}}(\mathrm{s}{Y}_{\mathrm{up}}-\mathrm{i}{Y}_{\mathrm{up}})dX,$$

(A8b)

and

$$X_0=\mathrm{s}{Y}_{\mathrm{down}}\cap \mathrm{i}{Y}_{\mathrm{up}}.$$

(A8c)

Figure A2. Area of the projection of the whole vertical wall on the plane of the sky $XY$, defined by the intersection of the ellipses $Y_{\mathrm{up}}$ and $Y_{\mathrm{down}}$.

Figure A2. Area of the projection of the whole vertical wall on the plane of the sky $XY$, defined by the intersection of the ellipses $Y_{\mathrm{up}}$ and $Y_{\mathrm{down}}$.

Since

$$\mathrm{s}{Y}_{\mathrm{up}}=R_{\mathrm{wall}}cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}}}\right)}^2}+H_{\mathrm{wall}}sin\left(i\right),$$

$$\mathrm{i}{Y}_{\mathrm{up}}=-R_{\mathrm{wall}}cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}}}\right)}^2}+H_{\mathrm{wall}}sin\left(i\right),$$

and

$$\mathrm{s}{Y}_{\mathrm{down}}=R_{\mathrm{wall}}cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}}}\right)}^2}-H_{\mathrm{wall}}sin\left(i\right),$$

it follows

$$A_1=2{\int }_0^{X_0}2H_{\mathrm{wall}}sin\left(i\right)dX=4X_0{H}_{\mathrm{wall}}sin\left(i\right),$$

$$A_2=2R_{\mathrm{wall}}cos\left(i\right){\int }_{X_0}^{R\mathrm{wall}}\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}}}\right)}^2}dX$$

and

$$X_0=R_{\mathrm{wall}}\sqrt{1-{\left[{\displaystyle \frac{H_{\mathrm{wall}}}{R_{\mathrm{wall}}}}tan\left(i\right)\right]}^2}.$$

Using Equation (A6), it leads

$$A_2=R_{\mathrm{wall}}cos\left(i\right)\left[\frac{1}{2}\pi R_{\mathrm{wall}}-X_0\sqrt{1-{\left({\displaystyle \frac{X_0}{R_{\mathrm{wall}}}}\right)}^2}-R_{\mathrm{wall}}{sin}^{-1}\left({\displaystyle \frac{X_0}{R_{\mathrm{wall}}}}\right)\right].$$

Hence, the area of the projected vertical wall is

$$\begin{array}{cc}\hfill A=& 4X_0{H}_{\mathrm{wall}}sin\left(i\right)\hfill \\ \hfill +& R_{\mathrm{wall}}cos\left(i\right)\left[\pi R_{\mathrm{wall}}-2X_0\sqrt{1-{\left({\displaystyle \frac{X_0}{R_{\mathrm{wall}}}}\right)}^2}-2R_{\mathrm{wall}}{sin}^{-1}\left({\displaystyle \frac{X_0}{R_{\mathrm{wall}}}}\right)\right].\hfill \end{array}$$

(A11)

Appendix A.4. Theoretical Area of the Projected Vertical Wall: The Wall with Shadow

To calculate the area of the projection of a vertical wall considering the umbra cast by the inner disk, I subtract to the area generated by the whole wall [see Equation (A11)], the area generated by the shadow. The shadow can be assumed as a vertical wall with $H_{\mathrm{shw}}$ in height, measured from the mid–plane.

Then the area of the projection of the shadow on the plane of the sky is

$$A_{\mathrm{shw}}=A_1+2A_2$$

(A12)

as seen in Figure A3. Where, according to Definition A1,

$$A_1={\int }_{-X_1}^{X_1}(\mathrm{s}{Y}_{\mathrm{up}}^{\prime }-\mathrm{s}{Y}_{\mathrm{down}}^{\prime })dX,$$

(A13a)

$$A_2={\int }_{X_1}^{X_2}(\mathrm{s}{Y}_{\mathrm{up}}^{\prime }-\mathrm{i}{Y}_{\mathrm{up}})dX,$$

(A13b)

and

$$X_1=\mathrm{i}{Y}_{\mathrm{up}}\cap \mathrm{s}{Y}_{\mathrm{down}}^{\prime }.$$

(A13c)

$$X_2=\mathrm{i}{Y}_{\mathrm{up}}\cap \mathrm{s}{Y}_{\mathrm{up}}^{\prime }.$$

(A13d)

with

$$0<X_0<X_1<X_2<R_{\mathrm{wall}}$$

Since

$$\mathrm{s}{Y}_{\mathrm{up}}^{\prime }=R_{\mathrm{wall}}cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}}}\right)}^2}+H_{\mathrm{shw}}sin\left(i\right),$$

$$\mathrm{s}{Y}_{\mathrm{down}}^{\prime }=R_{\mathrm{wall}}cos\left(i\right)\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}}}\right)}^2}-H_{\mathrm{shw}}sin\left(i\right),$$

if follows

$$A_1=2{\int }_0^{X_1}2H_{\mathrm{shw}}sin\left(i\right)dX=4X_1{H}_{\mathrm{shw}}sin\left(i\right),$$

$$A_2=2R_{\mathrm{wall}}cos\left(i\right){\int }_0^{X_1}\sqrt{1-{\left({\displaystyle \frac{X}{R_{\mathrm{wall}}}}\right)}^2}dX+{\int }_{X_1}^{X_2}(H_{\mathrm{shw}}-H_{\mathrm{wall}})sin\left(i\right)dX,$$

$$X_1=R_{\mathrm{wall}}\sqrt{1-{\left[{\displaystyle \frac{H_{\mathrm{wall}}+H_{\mathrm{shw}}}{2R_{\mathrm{wall}}}}tan\left(i\right)\right]}^2}$$

and

$$X_2=R_{\mathrm{wall}}\sqrt{1-{\left[{\displaystyle \frac{H_{\mathrm{wall}}-H_{\mathrm{shw}}}{2R_{\mathrm{wall}}}}tan\left(i\right)\right]}^2}$$

Figure A3. Area of the projection of the shadow cast by the inner disk on the plane of the sky $XY$, defined by the intersection of the ellipses $Y_{\mathrm{up}}^{\prime }$, $Y_{\mathrm{down}}^{\prime }$ and $Y_{\mathrm{up}}$ (in gray). The region in blue is the area of the visible wall.

Figure A3. Area of the projection of the shadow cast by the inner disk on the plane of the sky $XY$, defined by the intersection of the ellipses $Y_{\mathrm{up}}^{\prime }$, $Y_{\mathrm{down}}^{\prime }$ and $Y_{\mathrm{up}}$ (in gray). The region in blue is the area of the visible wall.

Using Equation (A6), it leads

$$\begin{array}{cc}\hfill A_2=& (H_{\mathrm{shw}}-H_{\mathrm{wall}})(X_2-X_1)sin\left(i\right)\hfill \\ \hfill +& R_{\mathrm{wall}}cos\left(i\right)\left[X_2\sqrt{1-{\left({\displaystyle \frac{X_2}{R_{\mathrm{wall}}}}\right)}^2}+R_{\mathrm{wall}}{sin}^{-1}\left({\displaystyle \frac{X_2}{R_{\mathrm{wall}}}}\right)\right.\hfill \\ \hfill -& \left.X_1\sqrt{1-{\left({\displaystyle \frac{X_1}{R_{\mathrm{wall}}}}\right)}^2}-R_{\mathrm{wall}}{sin}^{-1}\left({\displaystyle \frac{X_1}{R_{\mathrm{wall}}}}\right)\right].\hfill \end{array}$$

Hence, the area of the shadow, projected on the plane of the sky, is

$$\begin{array}{cc}\hfill A_{\mathrm{shw}}=& 4X_1{H}_{\mathrm{shw}}sin\left(i\right)+2(H_{\mathrm{shw}}-H_{\mathrm{wall}})(X_2-X_1)sin\left(i\right)\hfill \\ \hfill +& 2R_{\mathrm{wall}}cos\left(i\right)\left[X_2\sqrt{1-{\left({\displaystyle \frac{X_2}{R_{\mathrm{wall}}}}\right)}^2}+R_{\mathrm{wall}}{sin}^{-1}\left({\displaystyle \frac{X_2}{R_{\mathrm{wall}}}}\right)\right.\hfill \\ \hfill -& \left.X_1\sqrt{1-{\left({\displaystyle \frac{X_1}{R_{\mathrm{wall}}}}\right)}^2}-R_{\mathrm{wall}}{sin}^{-1}\left({\displaystyle \frac{X_1}{R_{\mathrm{wall}}}}\right)\right].\hfill \end{array}$$

(A16)