Light radiation pressure upon an optically orthotropic surface

TL;DR

This paper develops a model to calculate light radiation pressure on anisotropic, orthotropic surfaces, such as wrinkled solar sails, using modified Maxwell's equations and Bayesian analysis to improve accuracy over isotropic models.

Contribution

It introduces an orthotropic model for optical parameters and derives force equations considering emission, absorption, and reflection, enhancing modeling of complex surfaces.

Findings

Orthotropic model can outperform isotropic models in certain cases.

Bayesian analysis confirms the model's improved accuracy for wrinkled surfaces.

The model is applicable to flat solar sails with wrinkles.

Abstract

In this paper, we discuss the problem of determination of light radiation pressure force upon an anisotropic surface. The optical parameters of such a surface are considered to have major and minor axes, so the model is called an orthotropic model. We derive the equations for force components from emission, absorption, and reflection, utilizing a modified Maxwell's specular-diffuse model. The proposed model can be used to model a flat solar sail with wrinkles. By performing Bayesian analysis for example of a wrinkled surface, we show that there are cases in which an orthotropic model of the optical parameters of a surface may be more accurate than an isotropic model.

Tables

Table 1. Table 1. Summary of the posterior distributions for different simulations. B ​ F 𝐵 𝐹 BF – Bayesian Factor – posterior probability of the isotropic model divided by posterior probability of the orthotropic model. Parameters – median values of parameters for the most probable model with corresponding 95% Highest Density Interval (HDI) limits. Deviation – calculated mean-square deviation between Monte Carlo simulated data and calculated data based on median predicted values of the parameters.

Case	$BF$	Parameters	Deviation
Diffuse plate,	$\approx$0	${\rho }_1={0.443}_{-0.097}^{+0.148}$
${\rho }_0=1,s_0=0,B_0=2/3$		${\rho }_2={0.521}_{-0.130}^{+0.161}$
		$s_1={0.701}_{-0.201}^{+0.299}$
		$s_2={0.598}_{-0.209}^{+0.360}$
		$k={0.818}_{-0.073}^{+0.099}$
		${\sigma }_F={0.107}_{-0.009}^{+0.011}$	0.047
Specular plate,	297	$\rho ={0.999}_{-0.001}^{+0.001}$
${\rho }_0=1,s_0=1$		$s={0.999}_{-0.003}^{+0.001}$
		$B={0.663}_{-0.577}^{+0.336}$
		${\sigma }_F={0.0148}_{-0.0015}^{+0.0015}$	0.0108
		(e.g. Fig. 5)
Specular-Diffuse plate,	99	$\rho ={0.909}_{-0.085}^{+0.090}$
${\rho }_0=1,s_0=0.5,B_0=2/3$		$s={0.511}_{-0.058}^{+0.062}$
		$B={0.790}_{-0.131}^{+0.203}$
		${\sigma }_F={0.0742}_{-0.0072}^{+0.0082}$	0.0379

Table 2. Table 2. Summary of posterior distributions for different simulations for surface with waves with different optical parameters ( case column). Isotropic – summary of MCMC approximation considering an isotropic model of optical parameters, value is a median, limits are 95% HDI. Orthotropic – summary of MCMC approximation considering an orthotropic model of optical parameters, value is a median, limits are 95% HDI. Deviation – calculated mean-square deviation between Monte Carlo simulated data and calculated data based on median predicted values of the parameters.

Case	Isotropic	Orthotropic
${\rho }_0=1,s_0=1$	$\rho ={0.870}_{-0.105}^{+0.130}$	${\rho }_1={0.749}_{-0.156}^{+0.233}$
	$s={0.516}_{-0.136}^{+0.144}$	${\rho }_2={0.729}_{-0.174}^{+0.219}$
	$B={0.770}_{-0.200}^{+0.230}$	$s_1={0.803}_{-0.256}^{+0.191}$
	${\sigma }_F={0.281}_{-0.028}^{+0.030}$	$s_2={0.822}_{-0.297}^{+0.178}$
		$k={0.186}_{-0.085}^{+0.103}$
		${\sigma }_F={0.274}_{-0.024}^{+0.025}$
Deviation:	$0.108$	$0.093$
${\rho }_0=1,s_0=0.5,B_0=2/3$	$\rho ={0.574}_{-0.146}^{+0.325}$	${\rho }_1={0.291}_{-0.195}^{+0.323}$
	$s={0.010}_{-0.010}^{-0.034}$	${\rho }_2={0.476}_{-0.191}^{+0.324}$
	$B={0.776}_{-0.281}^{+0.224}$	$s_1={0.179}_{-0.179}^{+0.452}$
	${\sigma }_F={0.163}_{-0.017}^{+0.017}$	$s_2={0.222}_{-0.154}^{+0.394}$
		$k={1.870}_{-0.850}^{+3.850}$
		${\sigma }_F={0.122}_{-0.011}^{+0.016}$
Deviation:	$0.091$	$0.047$
${\rho }_0=1,s_0=0,B_0=2/3$	$\rho ={0.611}_{-0.196}^{+0.327}$	${\rho }_1={0.104}_{-0.023}^{+0.034}$
	$s={0.005}_{-0.005}^{+0.018}$	${\rho }_2={0.369}_{-0.085}^{+0.083}$
	$B={0.707}_{-0.240}^{+0.293}$	$s_1={0.677}_{-0.240}^{+0.308}$
	${\sigma }_F={0.190}_{-0.018}^{+0.021}$	$s_2={0.331}_{-0.136}^{+0.197}$
		$k={2.440}_{-0.510}^{+0.790}$
		${\sigma }_F={0.097}_{-0.009}^{+0.009}$
(e.g. Fig. 7)
Deviation:	$0.101$	$0.039$
${\rho }_0=0.5,s_0=1$	$\rho ={0.314}_{-0.169}^{+0.488}$	${\rho }_1={0.047}_{-0.047}^{+0.103}$
	$s={0.063}_{-0.063}^{+0.160}$	${\rho }_2={0.211}_{-0.123}^{+0.082}$
	$B={0.479}_{-0.337}^{+0.453}$	$s_1={0.094}_{-0.94}^{+0.462}$
	${\sigma }_F={0.139}_{-0.013}^{+0.015}$	$s_2={0.296}_{-0.207}^{+0.406}$
		$k={2.330}_{-1.130}^{+1.680}$
		${\sigma }_F={0.123}_{-0.012}^{+0.013}$
Deviation:	$0.068$	$0.040$
${\rho }_0=0.5,s_0=0.5,B_0=2/3$	$\rho ={0.332}_{-0.175}^{+0.505}$	${\rho }_1={0.057}_{-0.057}^{+0.077}$
	$s={0.008}_{-0.008}^{+0.039}$	${\rho }_2={0.074}_{-0.031}^{+0.133}$
	$B={0.519}_{-0.334}^{+0.456}$	$s_1={0.218}_{-0.206}^{+0.539}$
	${\sigma }_F={0.120}_{-0.012}^{+0.013}$	$s_2={0.448}_{-0.346}^{+0.456}$
		$k={4.450}_{-1.120}^{+2.130}$
		${\sigma }_F={0.081}_{-0.008}^{+0.007}$
Deviation:	$0.074$	$0.037$
${\rho }_0=0.5,s_0=0,B_0=2/3$	$\rho ={0.400}_{-0.192}^{+0.482}$	${\rho }_1={0.041}_{-0.013}^{+0.025}$
	$s={0.006}_{-0.006}^{+0.024}$	${\rho }_2={0.121}_{-0.036}^{+0.038}$
	$B={0.580}_{-0.346}^{+0.387}$	$s_1={0.806}_{-0.264}^{+0.194}$
	${\sigma }_F={0.137}_{-0.013}^{+0.015}$	$s_2={0.524}_{-0.206}^{+0.146}$
		$k={3.480}_{-1.170}^{+1.230}$
		${\sigma }_F={0.073}_{-0.006}^{+0.008}$
Deviation:	$0.083$	$0.031$

Equations

$${ϵ}^{\prime }(\beta ,\theta ,T)=\frac{\pi \underset{0}{\overset{\infty }{\int }}{ϵ}_{\lambda }^{\prime }{i}_{\lambda b}^{\prime }d\lambda }{\sigma T^4},$$

$$i_{\lambda b}^{\prime }(\lambda ,T)=\frac{2h{c}^2}{{\lambda }^5\left(e^{\frac{hc}{\lambda kT}}-1\right)},$$

$$\hat{\mathbf{r}}(\beta ,\theta )=\sin \beta \cos \theta {\hat{\mathbf{e}}}_1+\sin \beta \sin \theta {\hat{\mathbf{e}}}_2+\cos \beta {\hat{\mathbf{e}}}_3,$$

$$d{\mathbf{F}}^{Sr}(\beta ,\theta ,T)=-\frac{{ϵ}^{\prime }\sigma T^4}{c}\hat{\mathbf{r}}\cos \beta dA.$$

$$d{F}_i^S(T)=\frac{1}{2\pi }\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\frac{\pi }{2}}{\int }}d{\mathbf{F}}^{Sr}\cdot {\hat{\mathbf{e}}}_{i}d\beta d\theta dA.$$

$${ϵ}^{\prime }=ϵf_{ϵ}(\beta ),$$

$$d{F}_1^S=-\frac{\sigma T^4}{2\pi c}ϵ\underset{0}{\overset{\frac{\pi }{2}}{\int }}{f}_{ϵ}(\beta )\sin \beta \cos \beta \left(\underset{0}{\overset{2\pi }{\int }}\cos \theta d\theta \right)d\beta dA=0;$$

$$d{F}_2^S=-\frac{\sigma T^4}{2\pi c}ϵ\underset{0}{\overset{\frac{\pi }{2}}{\int }}{f}_{ϵ}(\beta )\sin \beta \cos \beta \left(\underset{0}{\overset{2\pi }{\int }}\sin \theta d\theta \right)d\beta dA=0;$$

$$d{F}_3^S=-\frac{\sigma T^4}{2\pi c}ϵ\underset{0}{\overset{\frac{\pi }{2}}{\int }}{f}_{ϵ}(\beta ){\cos }^2\beta \left(\underset{0}{\overset{2\pi }{\int }}d\theta \right)d\beta dA=-\frac{\sigma T^4}{c}ϵ\underset{0}{\overset{\frac{\pi }{2}}{\int }}{f}_{ϵ}(\beta ){\cos }^2\beta d\beta dA.$$

$$B=\underset{0}{\overset{\frac{\pi }{2}}{\int }}{f}_{ϵ}(\beta ){\cos }^2\beta d\beta .$$

$$B=\underset{0}{\overset{\frac{\pi }{2}}{\int }}{\cos }^3\beta d\beta =\frac{2}{3}.$$

$$d{\mathbf{F}}^S=-\frac{ϵB\sigma T^4}{c}\hat{\mathbf{n}}.$$

$${ϵ}_m^{\prime }(\beta ,\theta )=({ϵ}_1{\cos }^2(\theta -{\theta }_m)+{ϵ}_2{\sin }^2(\theta -{\theta }_m))\cos \beta .$$

$$d{F}_1^{Sm}=-\frac{\sigma T^4}{2\pi c}\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\frac{\pi }{2}}{\int }}{ϵ}_m^{\prime }(\beta ,\theta )\cos \beta \sin \beta \cos \theta d\beta d\theta dA=0;$$

$$d{F}_2^{Sm}=-\frac{\sigma T^4}{2\pi c}\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\frac{\pi }{2}}{\int }}{ϵ}_m^{\prime }(\beta ,\theta )\cos \beta \sin \beta \sin \theta d\beta d\theta dA=0;$$

$$d{F}_3^{Sm}=-\frac{\sigma T^4}{2\pi c}\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\frac{\pi }{2}}{\int }}{ϵ}_m^{\prime }(\beta ,\theta ){\cos }^2\beta d\beta d\theta dA=-\frac{\sigma T^4({ϵ}_1+{ϵ}_2)}{3c}dA.$$

$$B_m=\frac{{ϵ}_1+{ϵ}_2}{3{ϵ}_1}.$$

$$d{\mathbf{F}}^S=-\frac{{ϵ}_1{B}_m\sigma T^4}{c}\hat{\mathbf{n}}$$

$$i^{\prime A}(\beta ,\theta )=\underset{0}{\overset{\infty }{\int }}{i}_{\lambda }^{\prime A}d\lambda .$$

$$d{\mathbf{F}}^{Ar}(\beta ,\theta )=-\frac{i^{\prime A}}{c}\hat{\mathbf{r}}\cos \beta dA.$$

$$d{F}_i^A=\frac{1}{2\pi }\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\frac{\pi }{2}}{\int }}d{\mathbf{F}}^{Ar}\cdot {\hat{\mathbf{e}}}_{i}d\beta d\theta .$$

$$d{F}_i^A=-hc\underset{0}{\overset{\infty }{\int }}\underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\frac{\pi }{2}}{\int }}\frac{\hat{\mathbf{r}}\cdot {\hat{\mathbf{e}}}_i{i}_{\lambda }^{\prime A}}{2\pi h{c}^2}\cos \beta d\beta d\theta d\lambda dA.$$

$$i_{\lambda }^{\prime A}=2\pi q(\lambda )\cos \beta \delta (\beta -{\beta }_0)\delta (\theta -{\theta }_0),$$

$$q_0=\underset{0}{\overset{\infty }{\int }}q(\lambda )d\lambda .$$

$$d{F}_1^A=-\frac{q_0}{c}\underset{0}{\overset{\frac{\pi }{2}}{\int }}\underset{0}{\overset{2\pi }{\int }}\sin \beta \cos \beta \cos \theta \delta (\beta -{\beta }_0)\delta (\theta -{\theta }_0)d\theta d\beta dA=$$

$$=-\frac{q_0}{c}\sin {\beta }_0\cos {\beta }_0\cos {\theta }_{0}dA;$$

$$d{F}_2^A=-\frac{q_0}{c}\underset{0}{\overset{\frac{\pi }{2}}{\int }}\underset{0}{\overset{2\pi }{\int }}\sin \beta \cos \beta \sin \theta \delta (\beta -{\beta }_0)\delta (\theta -{\theta }_0)d\theta d\beta dA=$$

$$=-\frac{q_0}{c}\sin {\beta }_0\cos {\beta }_0\sin {\theta }_{0}dA;$$

$$d{F}_3^A=-\frac{q_0}{c}\underset{0}{\overset{\frac{\pi }{2}}{\int }}\underset{0}{\overset{2\pi }{\int }}{\cos }^2\beta \delta (\beta -{\beta }_0)\delta (\theta -{\theta }_0)d\theta d\beta dA=-\frac{q_0}{c}{\cos }^2{\beta }_{0}dA.$$

$$\hat{\mathbf{s}}({\beta }_0,{\theta }_0)=-\sin {\beta }_0\cos {\theta }_0{\hat{\mathbf{e}}}_1-\sin {\beta }_0\sin {\theta }_0{\hat{\mathbf{e}}}_2-\cos {\beta }_0{\hat{\mathbf{e}}}_3,$$

$$d{\mathbf{F}}^A=-\frac{q_0}{c}(\hat{\mathbf{n}}\cdot \hat{\mathbf{s}})\hat{\mathbf{s}}dA.$$

$${\rho }^{\prime \prime }({\beta }^R,{\theta }^R,\beta ,\theta )=\frac{\underset{0}{\overset{\infty }{\int }}{\rho }_{\lambda }^{\prime \prime }{i}_{\lambda }^{\prime A}d\lambda }{i^{\prime A}},$$

$${\hat{\mathbf{r}}}^R(\beta ,\theta )=\sin {\beta }^R\cos {\theta }^R{\hat{\mathbf{e}}}_1+\sin {\beta }^R\sin {\theta }^R{\hat{\mathbf{e}}}_2+\cos {\beta }^R{\hat{\mathbf{e}}}_3.$$

$$I^{\prime }=\frac{1}{2\pi }\underset{0}{\overset{\frac{\pi }{2}}{\int }}\underset{0}{\overset{2\pi }{\int }}{\rho }^{\prime \prime }{i}^{\prime A}d\theta d\beta .$$

Full text

Nomenclature

$O^{\prime}x^{\prime}_{1}x^{\prime}_{2}x^{\prime}_{3}$ – global coordinate frame

$Ox_{1}x_{2}x_{3}$ – local coordinate frame

$\hat{\mathbf{e}}^{\prime}_{i},\ i=1,2,3$ – unit vectors of global coordinate frame

$\hat{\mathbf{e}}_{i},\ i=1,2,3$ – unit vectors of local coordinate frame

$dA$ – infinitesimal element of surface $A$

$\hat{\mathbf{n}}$ – normal to $dA$

$\hat{\mathbf{m}}$ – orientation vector for orthotropic model (in plane $Ox_{1}x_{2}$)

$\mathbf{r}$ – position of $dA$ in global frame

$\theta,\beta$ – direction angles in local frame

$\epsilon^{\prime}_{\lambda}$ – directional spectral emissivity

$\epsilon^{\prime}$ – directional integral emissivity

$\epsilon$ – emissivity (for isotropic case)

$B$ – Lambertian coefficient

$T$ – temperature of $dA$

$\epsilon_{1},\ \epsilon_{2},\ \theta_{m}$ – parameters of orthotropic model for emission

$B_{m}$ – modified Lambertian coefficient for orthotropic emission

$c$ – speed of light in vacuum

$i^{\prime A}_{\lambda}$ – directional spectral intensity of irradiation

$i^{\prime A}$ – directional integral intensity of irradiation

$q_{0}$ – integral intensity of light source

$\hat{\mathbf{s}}$ – vector from light source to $dA$

$\rho^{\prime\prime}_{\lambda}$ – bidirectional spectral reflectivity

$\rho^{\prime\prime}$ – bidirectional integral reflectivity

$I^{\prime}$ – hemispherical-directional light intensity

$s$ – specularity coefficient

$\rho$ – reflectivity (for isotropic model)

$\rho_{1},\ \rho_{2},\ \theta_{m}$ – parameters for orthotropic model for reflection

$B_{\rho}$ – modified Lambertian coefficient for orthotropic reflection

$d\mathbf{F}^{Sr}$ – fraction of emission pressure in arbitrary direction $\hat{\mathbf{r}}$

$d\mathbf{F}^{Ar}$ –fraction of absorption pressure in arbitrary direction $\hat{\mathbf{r}}$

$d\mathbf{F}^{Rr}$ – fraction of reflection pressure in arbitrary direction $\hat{\mathbf{r}}^{R}$

$d\mathbf{F}^{S}$ – light pressure from emission

$d\mathbf{F}^{A}$ – light pressure from absorption

$d\mathbf{F}^{R}$ – light pressure from reflection

$d\mathbf{F}$ – total light radiation pressure upon $dA$

Introduction

The theory of light radiation pressure upon space objects is well-established. For celestial bodies, this pressure creates the Yarkovsky acceleration due to uneven heating of their surface Vokrouhlický and Farinella (1998); Hartmann et al. (1999). There is also a Yarkovsky-O’Keefe–Radzievskii–Paddack (YORP) effect, in which an asteroid can spin-up from emission pressure because of its irregular shape, Vokrouhlický and Čapek (2002) up to the disintegration of a body Paddack (1969); Rubincam (2000).

For practical applications, the derivation of light radiation pressure force is necessary for the prediction of the dynamics of GNSS satellites Fliegel and Gallini (1996); Bar-Sever and Russ (1997); Springer et al. (1999); Bar-Sever and Kuang (2004); Rodriguez-Solano et al. (2012); Tan et al. (2016), for interplanetary stations Kubo-oka and Sengoku (1999); Vaughan et al. (2001); Turyshev et al. (2012) and other spacecraft Kinzel (2010).

For solar sail applications, there are many studies of light radiation pressure, including light pressure generalizations Forward (1989); McInnes (2004), and special cases – variable reflectance / transmittance coatings Kislov (2004), degradation effects Dachwald et al. (2005, 2006), joint analysis of aerodynamic and radiation forces on spacecraft Shmatov and Mordvinkin (2014), laser propulsion Forward (1984); Popova et al. (2016), transparent sails Swartzlander Jr (2017), etc. There are numerous studies of the astrodynamics of solar sails Farrés and Jorba (2016); Gachet et al. (2016); Lachut and Bennett (2016); Ono et al. (2016); Felicetti et al. (2017); German et al. (2017); Ma et al. (2017); Niccolai et al. (2017) etc.

In the space experiments Nanosail-D2 Alhorn et al. (2011), IKAROS Tsuda et al. (2011); Kawaguchi (2014), and LightSail Ridenoure et al. (2015) it was shown that any solar sail membrane has general curvature, both regular and semi-random (smoothness), and also small wrinkles.

The light pressure model on the curved solar sail was generalized by Rios-Reyes and Scheeres (2004, 2005); Rios-Reyes (2006); Rios-Reyes and Scheeres (2007); Scheeres (2007); McMahon and Scheeres (2010, 2014); McMahon and Scheeres (2015) and extended by Jing et al. (2012, 2014); Nerovny et al. (2017). This model is called the Generalized Sail Model (GSM).

In this paper, we will consider the optical anisotropy from the geometrical sources of this anisotropy. The main sources of this optical anisotropy are wrinkles on the solar sail membrane Wong and Pellegrino (2006a); Jenkins (2006). One special case of the effects of wrinkles on the solar sail efficiency was studied by Greschik (2014).

We will derive the equations for light radiation pressure by utilizing the well-established theory of light-matter interaction as in radiative heat transfer Howell et al. (2015), and after this, we will move to the vector representation of force.

We will consider the effects of emission, absorption, and reflection on light pressure because further phenomena such as transmission are supposed to be less influential on solar sails than the main effects Forward (1989). For each effect, we will consider both isotropic and anisotropic cases. For reflection, we will utilize Maxwell’s reflection model Howell et al. (2015), in which we assume that reflection has two components as a sum of diffuse and specular cases with corresponding specularity coefficient $s$. We will also consider the back reflection phenomenon for the orthotropic model.

1. Reference frames

Let us consider some surface $A$ in Euclidean space with origin $O^{\prime}$ and associated Cartesian coordinate system $O^{\prime}x^{\prime}_{1}x^{\prime}_{2}x^{\prime}_{3}$, Fig. 1. We will call this frame a global frame. Let us introduce $\hat{\mathbf{e}}^{\prime}_{i}$ – unit vectors for the global frame, $i=1{,}\ 2{,}\ 3$.

On this surface, it is possible to localize an infinitesimal surface element $dA$ for which we introduce a local Cartesian coordinate system $Ox_{1}x_{2}x_{3}$ with unit vectors $\hat{\mathbf{e}}_{i}$, Fig. 1. The origin of local frame $O$ is situated in the center of $dA$, and its normal $\hat{\mathbf{n}}$ is equal to $\hat{\mathbf{e}}_{3}$. $[T]$ is a transformation matrix from the local frame to the global frame. The orientation of $Ox_{1}$ and $Ox_{2}$ is arbitrary.

In the following equations, for any vector, e.g. $\hat{\mathbf{r}}$, we will use direction angles $(\beta,\theta)$ in the local frame as follows (Fig. 2):

• •

$\beta\in[0,\pi/2]$ – angle between vector $\hat{\mathbf{r}}$ and $+x_{3}$. We consider that the infinitesimal surface element is laying on the plane $Ox_{1}x_{2}$.

• •

$\theta\in[0,2\pi]$ – angle between axis $Ox_{1}$ and a projection of $\hat{\mathbf{r}}$ on the plane $Ox_{1}x_{2}$, counterclockwise around $Ox_{3}$.

Direction angles $(\beta,\theta)$ may have additional subscripts or superscripts.

2. Model

2.1. Thermal emission

Let $\epsilon^{\prime}_{\lambda}$ be a directional spectral emissivity, which depends on wavelength, temperature, and shows the difference of emission of $dA$ in direction $(\beta,\theta)$ as compared with black body emission in the same direction. One can write the equation of directional integral emissivity Howell et al. (2015):

[TABLE]

Where $\sigma$ – Stephan-Boltzmann constant, $i^{\prime}_{\lambda b}(\lambda,T)$ – spectral intensity of blackbody radiation which is represented by Planck’s law:

[TABLE]

Where $h$ – Planck’s constant, $c$ – light speed in vacuum, $k$ – Boltzmann constant.

Let us introduce an arbitrary unit vector $\hat{\mathbf{r}}$ in the local frame:

[TABLE]

Where $\beta\in[0;\pi/2]$ and $\theta\in[0;2\pi]$.

One can write the equation of the fraction of light radiation pressure in direction $\hat{\mathbf{r}}$:

[TABLE]

The superscript $S$ stands for emission (“self”). The relation $\cos\beta dA$ is the area of the infinitesimal element $dA$ under the angle $\beta$.

The equation for projection of light pressure force from emission in the direction $\hat{\mathbf{e}}_{i}$ on the area $dA$ can be written as follows:

[TABLE]

Isotropic case

Let us assume that the optical parameters are independent of direction within the surface i.e. they are isotropic. Function $\epsilon^{\prime}$ can be considered as an axis-symmetrical function around $Ox_{3}$,

[TABLE]

Where $\epsilon=\text{const}$ – emissivity of the material.

After rewriting of (1) in corresponding projections, using (2), we can obtain the following relations:

[TABLE]

Now we introduce the coefficient $B$:

[TABLE]

This value represents the composition of the axially symmetric radiation pattern of emission. It is often called a Lambertian coefficient Jing et al. (2012), since its most common value $B=2/3$ corresponds to a Lambertian diffuse surface. For diffuse surface $f_{\epsilon}(\beta)=\cos\beta$ (Lambert’s law),

[TABLE]

We can rewrite the equation for light pressure force from emission using vector notation:

[TABLE]

Eq. (3) is independent of the coordinate frame choice, since we can always represent $\hat{\mathbf{n}}$ in different coordinate systems using the corresponding transformation matrix $[T]$.

Orthotropic case

We will utilize the following model of optical parameters. There are two axes of optical parameters. In the direction that is defined by angle $\theta_{m}$ in the plane $Ox_{1}x_{2}$, the emissivity is equal to $\epsilon_{1}$. In the perpendicular direction, the emissivity is equal to $\epsilon_{2}$. We will assume that emissivity changes analytically according to some law (Fig. 3). So we can write the equation for directional integral emissivity as follows:

[TABLE]

For the isotropic case, assume $\epsilon_{1}=\epsilon_{2}=\epsilon$ for any $\theta_{m}$ $\epsilon^{\prime}_{m}(\beta,\theta)=\epsilon\cos\beta$.

Now we can write the projections of light pressure force from heat emission:

[TABLE]

Now we can introduce the modified Lambertian coefficient $B_{m}$ as follows:

[TABLE]

By introducing of the normal unit vector $\hat{\mathbf{n}}$, we can write the equation for emission pressure from an orthotropic surface:

[TABLE]

for isotropic case, where $\epsilon_{1}=\epsilon_{2}=\epsilon$, $B_{m}=2/3.$

2.2. Absorption

For light pressure from absorption, we will consider that all of the light flux linear momenta that fall onto $dA$ are transferred to that surface Forward (1989) so there is no difference between isotropic and anisotropic models. The anisotropy of absorption, however, should be considered in the case of calculation of solar sail temperature.

We can introduce the directional spectral intensity of irradiation $i^{\prime A}_{\lambda}(\lambda,\beta,\theta)$, which depends on the properties and spatial position of the light source relative to $dA$. We will also utilize the directional integral intensity of irradiation:

[TABLE]

Light pressure force fraction in direction $\hat{\mathbf{r}}$ can be written as follows:

[TABLE]

The superscript $A$ stands for the absorption.

One can write the equation for projection of light pressure force from absorbed radiation in the direction $\hat{\mathbf{e}}_{i}$:

[TABLE]

After substitution of Eq. (6) into Eq. (7), using Eq. (5), we can get:

[TABLE]

We will only consider the point light source with position $(\beta_{0},\theta_{0})$. The intensity of the falling light can be written as follows:

[TABLE]

Where $\delta(x)$ – Dirac’s delta function, $q(\lambda)$ – spectral intensity of light source for which the integral intensity is equal to $q_{0}$:

[TABLE]

The derivation for light pressure force from non-point light source can be found in McInnes and Brown (1990a, b).

After substitution of (9) and (10) into (8), we can get the projections of light pressure force for fully absorbed light:

[TABLE]

We can introduce the additional unit vector of light orientation $\hat{\mathbf{s}}$, which is pointing from light source to surface element $dA$:

[TABLE]

After this we can rewrite the equation for light pressure from absorbed light in simplified vector notation:

[TABLE]

2.3. Reflection

Now we will consider the light radiation pressure from reflected light. Let us introduce the bidirectional spectral reflectivity $\rho^{\prime\prime}_{\lambda}(\lambda,\beta^{R},\theta^{R},\beta,\theta)$, where $(\beta^{R},\theta^{R})$ is an orientation of reflected intensity. Considering the directional spectral intensity $i^{\prime A}_{\lambda}(\lambda,\beta,\theta)$, one can write the bidirectional integral reflectivity:

[TABLE]

Where $i^{\prime A}$ can be calculated using (5).

Let us introduce the arbitrary unit vector $\hat{\mathbf{r}}^{R}$, which represents the direction of reflected light flux:

[TABLE]

We can calculate the hemispherical-directional light intensity:

[TABLE]

The fraction of light pressure force from reflected light in direction $\hat{\mathbf{r}}^{R}$ will be as follows:

[TABLE]

and the total infinitesimal force from reflected light can be calculated for diffuse reflection as the follows:

[TABLE]

and for specular reflection,

[TABLE]

since all amount of specularly reflected flux will reflect in the same direction, as opposed to the diffuse case, in which we should take into account the projection of $dA$ in the direction perpendicular to $\hat{\mathbf{r}}^{R}$, i.e. $\cos\beta^{R}dA$.

2.3.1. Diffuse reflection

Isotropic case

For axially symmetric (diffuse) reflection the bidirectional reflectivity can be represented by the following equation:

[TABLE]

By use of Eq. (9) we can write the hemispherical-directional intensity:

[TABLE]

For infinitesimal force we can get the following relation:

[TABLE]

moreover, in coordinate projections we can get:

[TABLE]

By utilizing vector notation, after introducing the light source orientation vector $\hat{\mathbf{s}}$, we can obtain the equation of light pressure force for diffusely reflected light:

[TABLE]

Orthotropic case

For the orthotropic case the bidirectional reflectivity will be modeled as the following relation (Fig. 4):

[TABLE]

For specularity $s$ we will also utilize the similar orthotropic model:

[TABLE]

The hemispherical-directional intensity can be obtained by integration:

[TABLE]

After integration, we can get the following projections of force:

[TABLE]

Where

[TABLE]

For the isotropic case $\rho_{1}=\rho_{2}=\rho$, $s_{1}=s_{2}=s$, $B_{\rho}=B=2/3$.

The equation for light pressure force from diffusely reflected light from an orthotropic surface can be represented in the following vector form:

[TABLE]

2.3.2. Specular reflection

Isotropic case

For the specular case we can introduce the bidirectional reflectivity as follows:

[TABLE]

After integrating we can obtain the hemispherical-directional reflectivity:

[TABLE]

We can simply evaluate the following equations for projections of light pressure force:

[TABLE]

in vector notation

[TABLE]

Orthotropic case

For the orthotropic case the bidirectional reflectivity will be as follows:

[TABLE]

The projections of light pressure from specularly reflected light will be:

[TABLE]

We can introduce the unit vector of orientation of reflection axes in the local frame:

[TABLE]

After transformations, we can write the following equation for infinitesimal light pressure force from specularly reflected light from an optically orthotropic surface:

[TABLE]

For an orthotropic case of specular reflection, we also introduce the term from back reflection. The bidirectional reflectivity will be:

[TABLE]

Where $k$ – some empirical parameter.

The projections of light pressure from back reflected light will be:

[TABLE]

In vector notation:

[TABLE]

We included the back-reflection term according to numerical simulations of some wrinkled surface. The parameter $k$ of back reflection is a model constant, and it is not necessary that it have a strong physical background. In the Discussion section, we will obtain the value of $k$ simultaneously with the other optical parameters. We should note that for the isotropic model, $k$ should always be equal to 0. The absence of proper physical background for parameter $k$ is a major disadvantage of proposed model. However, it provides better coincidence with ray tracing results. See Discussion section for comparison between optically orthotropic model and ray tracing simulation for some wrinkled surface.

2.4. Equation of infinitesimal force and torque

Isotropic case

By substitution of terms from emission (3), absorption (11), diffuse (12) and specular reflection (14), we can obtain the total light pressure infinitesimal force of light pressure on element $dA$ for isotropic case:

[TABLE]

and, by transition from the local frame to the global frame using transformation matrix $[T]$, we can obtain the same equation for infinitesimal force in the global frame.

The obtained equation is similar to the solar sail light pressure force equation that is widely used Forward (1989); McInnes (2004); Rios-Reyes and Scheeres (2005); Jing et al. (2014); however we consider the emission only from the front surface. It is possible to introduce the emission pressure from the back side, i.e.

[TABLE]

However, it is not always necessary, because this can be a formulation of light pressure not only upon the solar sail but some structure with large internal volume (i.e. not a thin film). For solar sails, there is another approach, in which one can consider two opposite sides of the solar sail as two different surfaces very close to each other, introducing the visibility function (whether a particular side is illuminated or not). In this approach, it is not necessary to explicitly define which side is front and which is back, but the equation for light pressure force has much-complicated form Nerovny et al. (2017).

It is also possible to obtain the temperature $T$ of $dA$ using the thermal flux equilibrium equation and write the emission term without using $T$, e.g. Rios-Reyes (2006).

Orthotropic case

The equation for infinitesimal light pressure force can be summarized using equations (4), (11), (13), (15) and (16).

Light pressure torque

In both cases, isotropic and orthotropic, the infinitesimal light pressure torque can be obtained by the following cross product:

[TABLE]

Where $\mathbf{r}=(r^{\prime}_{1},r^{\prime}_{2},r^{\prime}_{3})^{T}$ – vector from the origin of the global frame $O^{\prime}x^{\prime}_{1}x^{\prime}_{2}x^{\prime}_{3}$ to the infinitesimal area $dA$.

3. Discussion

3.1. Bayesian analysis of parameters of the models

Bayesian analysis of parameters

In this section, we will utilize the so-called Bayesian analysis of model parameters Kruschke (2015); Kruschke and Liddell (2017). In this method, one has to follow several steps to determine which model is more probable than another and which parameters of models are more probable. These are the steps of Bayesian analysis Kruschke (2015):

1. (1)

Specify the experimental data, its scales, and identify which values should be predicted by the model and which data are predictors. 2. (2)

Define the descriptive analytical models of the data. 3. (3)

Specify the prior distribution for parameters. 4. (4)

Specify a likelihood function and use Bayesian inference to calculate the posterior probability distribution of parameter values. 5. (5)

Check that the posterior distribution approximates the data.

Data

For this analysis, we obtain the probability distribution of optical parameters of models for optically isotropic ($\rho,\ s,\ B$) and orthotropic surfaces ($\rho_{1},\ \rho_{2},\ s$). We do not consider light pressure force from emission. We do not divide values of force by speed of light $c$ and multiply them by flux $q_{0}$; thus, dimensions of force in all following tables and figures are $m^{2}$. To compare these models, we introduce the model index $m$. For the likelihood function, we also introduce the standard deviation $\sigma_{F}$.

Ray tracing

We wrote software for direct Monte Carlo simulation of ray tracing. This software represents the surface as a mesh of triangles. We utilize Maxwell’s model of specular-diffuse reflection. The diffuse reflection is Lambertian. This software does not consider spectral parameters of a light source. Calculation of light pressure force was done similar to works Nerovny et al. (2017); Ziebart (2004). We tested this software on several simple geometries which have analytical models of light radiation pressure: specular solar sail, Lambertian diffuse solar sail, specular and diffuse sphere, cylinder and cone. This software is freely available on GitHub Nerovny and Grigorjev (2017) under GPLv3 license.

Models

We use two models: a model of an optically isotropic surface and a model of an optically orthotropic surface, $\hat{\mathbf{m}}=(1,0,0)^{T}$ in local frame. Both models should predict vectors of light radiation pressure $\mathbf{F}_{\text{I}}$ and $\mathbf{F}_{\text{O}}$ consequently, for different orientations of a light source, and these vectors should be as close as possible to data $\mathbf{F}$, which is obtained by direct ray tracing.

Prior distribution

Prior distribution for the isotropic model: $\rho,\ s,\ B$ is uniform from 0 to 1, $\sigma_{F}$ is exponential with parameter equal to 1 Plummer (2013). Prior distribution for the orthotropic model: $\rho_{1},\ \rho_{2},\ s_{1},\ s_{2}$ is uniform from 0 to 1, $\sigma_{F}$ and $k$ are exponential with parameter equal to 1. Prior distribution for model index $m$: generally, for the isotropic model $P(m=1)=0.5$, for the orthotropic model $P(m=2)=0.5$. In subsequent analysis, we use different values for prior probabilities of a model index.

Likelihood function

The likelihood probability distribution function for this model is the following formula:

[TABLE]

Where: $C_{c}$ – an integer large enough to ensure that the likelihood is less than 1 (required by JAGS Plummer (2013)); $\mathbf{F}$ – data from direct Monte Carlo simulation (resultant light pressure force); $\mathbf{F}_{\text{I}}$ – predicted vector of light pressure force for given set of parameters considering the isotropic case, Eq. (17); $\mathbf{F}_{\text{O}}$ – predicted vector of light pressure force for given set of parameters considering orthotropic case, Eq. (11), (13) and (15); $\sigma_{F}$ – standard deviation.

Posterior distribution

In our analysis, we use the R code combined with the JAGS toolkit for Markov Chain Monte Carlo (MCMC) approximation of posterior probability Plummer (2013). The structure of this code is similar to examples in Kruschke (2015).

Check of prediction

For median values of parameters, we use the mean-square deviation between light pressure force predicted by the corresponding model and light pressure force from ray tracing.

3.1.1. Bayesian analysis for flat surface

Data

In this case, we analyzed two different models of light radiation pressure upon a flat surface. The surface is a square area in $O^{\prime}x^{\prime}_{1}x^{\prime}_{2}$ plane of the global coordinate system. The length of each side is equal to 1m, all sides are parallel to the corresponding axes of the coordinate system, $x^{\prime}_{1},\ x^{\prime}_{2}\in[-0.5,0.5],\ x^{\prime}_{3}=0$. Optical parameters of this surface are uniform and can be specular, diffuse or specular-diffuse with corresponding specularity coefficient: $\rho_{0}=1$ – reflectivity of the surface, $s_{0}\in\{0,0.5,1\}$ – specularity of the surface, $B_{0}=2/3$ – Lambertian coefficient of the front side, 3 combinations total. We calculated the light radiation pressure vector for all combinations of orientation angles of a light source from these limits:

• •

$\beta\in\left\{9\degree,18\degree,\dots,72\degree,81\degree\right\}$;

• •

$\theta\in\left\{0\degree,18\degree,\dots,342\degree,360\degree\right\}$,

189 combinations total.

The number of rays was 100000 in each simulation.

Posterior distribution

Tab. LABEL:tab:BF represent the summary of results of MCMC approximation of posterior distribution for different cases. In this table, the Bayesian Factor ($BF$) is the relative posterior probability of a number of a model which is most probable for given data. We should note that we used median values instead of the mode of distribution because the speed of convergence for the mode is much lower a than for median.

We performed the analysis for different prior probabilities of a models, e.g. for the specular case we set $P(m=1)=0.01$ and $P(m=2)=0.99$, the MCMC calculated value was $BF^{\prime}$, and, considering that both prior probabilities should equal to 0.5, according to Bayes rule, the final Bayesian Factor was Kruschke (2015):

[TABLE]

For the diffuse case, the MCMC approximation of posterior distribution showed that the optically orthotropic model of optical parameters is more probable than the isotropic, however, the median values for optical parameters represent the isotropic model very close. This problem arises when one tries to compare two models, one of which can represent another Kruschke (2015). In the subsequent analysis, we will not compare these models using the Bayesian Factor. Instead, we will compare these two models by calculating the mean-square deviation between actual data from direct Monte Carlo simulation of ray tracing and predicted data, which will be calculated using isotropic or orthotropic models and the corresponding median values of optical parameters.

3.1.2. Bayesian analysis for curved surface

Data

In this case, we analyzed two different models of light radiation pressure upon a surface with waves (wrinkles). As it was shown by different authors, both analytically Epstein (2003); Wong and Pellegrino (2006b), numerically Wong and Pellegrino (2006c); Wang et al. (2009); Xiao et al. (2011) and experimentally Wong and Pellegrino (2006a), the surface of wrinkled thin membrane can be represented according to the sine or cosine laws.

The magnitude of wrinkles can be calculated by the following formula Epstein (2003):

[TABLE]

Where: $\delta=\delta(x^{\prime}_{1},x^{\prime}_{2})$ – smooth field of amplitude of wrinkles over the wrinkled domain $\mathcal{D}$; $\lambda=\lambda(x^{\prime}_{1},x^{\prime}_{2})$ – smooth field of wavelength of wrinkles over the wrinkled domain $\mathcal{D}$.

We used the following parameters of wrinkled surface: $x^{\prime}_{1},\ x^{\prime}_{2}\in[-0.5,0.5]$, $x^{\prime}_{3}=\delta(\cos 8\pi x^{\prime}_{2}-1)$. There are 4 waves, wavelength $\lambda=const=0.25m$, amplitude $\delta=const=0.2m$, the global frame is shifted along $+x^{\prime}_{3}$ by $\delta/2$. This surface is shown in Fig. 6. Optical parameters of this surface are uniform and can be specular, diffuse or specular-diffuse with corresponding specularity coefficient: $\rho_{0}\in\{0.5,1\}$ – reflectivity of surface, $s_{0}\in\{0,0.5,1\}$ – specularity of surface, $B_{0}=2/3$ – Lambertian coefficient of the front side, 6 combinations total. We calculated the light radiation pressure vector for all combinations of orientation angles of a light source from these limits:

• •

$\beta\in\left\{9\degree,18\degree,\dots,72\degree,81\degree\right\}$;

• •

$\theta\in\left\{0\degree,18\degree,\dots,342\degree,360\degree\right\}$,

189 combinations total.

The number of rays was 100000 in each simulation. The ray tracing results are available in the Mendeley Data Nerovny (2017).

Posterior distribution

Tab. LABEL:tab:cosine represents the summary of the MCMC approximation of parameters of models for different values of $\rho_{0}$ and $s_{0}$. There are several cases for which the mean-square deviation of an orthotropic surface is less than for an isotropic surface; i.e., the orthotropic model is more accurate. This fact confirms that the model of an optically orthotropic surface proposed in this paper can be suitable for practical applications, providing a more accurate representation than an isotropic model, at least if not considering thermal emission.

3.2. Conformity with Generalized Sail Model

The other difficulty with a model of an optically orthotropic surface is the fact that it is not well fitted with the Generalized Sail Model Rios-Reyes and Scheeres (2005), since it uses an additional vector of an orientation of optical axes $\hat{\mathbf{m}}$, as soon as there is no such vector in the original GSM. One of the possible ways of dealing with this may be the expansion of terms with $\hat{\mathbf{m}}$ into some power series. Further analytical separation of $\hat{\mathbf{m}}$ from $\hat{\mathbf{s}}$ may be accomplished in the same way as it was done for $\hat{\mathbf{n}}$ and $\hat{\mathbf{s}}$ in Nerovny et al. (2017). Thus conformity with GSM needs to be investigated.

For practical applications without the GSM, it is possible to utilize the orthotropic model for better accuracy after derivation of model parameters for given surface.

4. Acknowledgements

The authors would like to thank assistant Dmitry A. Goncharov from the department “Theoretical Mechanics” of BMSTU and Evgeny S. Golubev from Astro Space Center of the P.N. Lebedev Physical Institute of the Russian Academy of Sciences for their valuable advice and discussions. The authors also would like to thank Mark A. Bowman and Eugene Chebezov from the Flight Dynamics Division, NASA Johnson Space Center, for their help in proofreading of the article.

Contributions: N.A. Nerovny – model and analysis, development of srp2 software, I.E. Lapina – Bayesian analysis, A.S. Grigorjev – srp2 software.

This work was done during the development of BMSTU-Sail Space Experiment Rachkin et al. (2011); bms (2014).