Consider the following stochastic dynamical system. 


Fix $a > 0$, $b > 0$, $c>0$ and $v > 0$, and let $\mathbf{r}(t)=(x(t),y(t),z(t))$ be the position at time $t$ of a point which moves in the parallelepiped $R=\{ (x,y,z) \in \mathbb{R}^3: 0 \leq x \leq a, 0 \leq y \leq b, 0 \leq z \leq c \}$ with velocity of constant magnitude $v$ a according to the following rules (we shall consider the usual spherical coordinates $(\rho, \theta, \phi)$ in $\mathbb{R}^3$, where $\theta$ is the polar angel with respect to the $z$-axis and $\phi$ the azimuthal angle):


(i) in the interior of $R$ the point is subject to no force, so that is moves with constant velocity $(v\_x(t),v\_y(t),v\_z(t))=\mathbf{v}(t)=\frac{d\mathbf{r}}{dt}(t)$;


(ii) when the point reaches the face of $R$ corresponding to $z=0$ it is reflected diffusely, that is its velocity after the collision has always magnitude $v$, and the versor $\mathbf{v} / v$ has polar coordinates $(1,\theta, \phi)$ with $(\theta, \phi)$ distributed over $[0,\pi/2] \times [0,2\pi)$ with probability density law (Knudsen law) $f(\theta, \phi)= \frac{1}{\pi} \cos \theta d \Omega= \frac{1}{\pi} \cos \theta \sin \theta d \theta d \phi$;


(iii) if the point reaches any other side of $R$, then it is reflected elastically: e.g. in the case of the face corresponding to $z=c$ the $x$ and $y$ components of $\mathbf{v}$ are preserved, while the $z$-component of $\mathbf{v}$ changes sign;


(iv) finally, if the point reaches an edge (or a vertex) of $R$, then the reflection must be thought as the combination of reflections of the two (respectively three) faces involved; so if the point hits the edge $\{ (x,0,0): 0 < x < a \}$ it is reflected diffusely, that is its velocity after the collision has always magnitude $v$, and the versor $\mathbf{v} / v$ has polar coordinates $(1,\theta, \phi)$ with $(\theta, \phi)$ distributed over $[0,\pi/2] \times [0, pi]$ with probability density law $g(\theta, \phi)= \frac{2}{\pi} \cos \theta d \Omega= \frac{2}{\pi} \cos \theta \sin \theta d \theta d \phi$; if instead e.g. the point hits the edge $\{ (x,0,a): 0 < x < a \}$ then is reflected elastically, that is the $x$ component is preserved, while the $y$ and $z$ components change sign.


For any time $T> 0$, let $N(T)$ be the number of times the point hits the face $z=c$ (included its edges and vertices) before or at time $T$. Let $(1,\theta\_j,\phi\_j)$ be the polar coordinates of the versor $\mathbf{v} / v$ just before the $j$-th collision and form the random sum
\begin{equation}
\sum\_{j=1}^{N(T)} \cos \theta\_j.
\end{equation}
I would like to prove that for any initial conditions $(\mathbf{r}(0), \mathbf{v}(0))$, with $v\_{z}(0) \neq 0$, the following equation holds
\begin{equation}
\lim\_{T \rightarrow \infty} \frac{1}{T} E \left[ \sum\_{j=1}^{N(T)} \cos \theta\_j \right] = \frac{v}{6 c}.
\end{equation}


NOTE1. This problem has been suggested to me by an interesting proof of a physical law called Wien's displacement law given by Richtmyer, Kennard and Cooper in their book "Introduction to Modern Physics", Sixth Edition, Appendix of Chapter 5 (see in particular p. 145, and see also ter Haar and Wergeland, Elements of Thermodynamics, $\S 5.3.3$). I formulated a first model of dynamic stochastic billiard in my post [Stochastic Dynamic Billiard](https://mathoverflow.net/questions/282145/a-stochastic-dynamical-billiard). This model was too complicated to be carefully analyzed, so I simplified it to get a much simpler version of it in a second post [A Simple Stochastic Dynamic Billiard](https://mathoverflow.net/questions/285729/a-simple-stochastic-dynamic-billiard). The answer by Iosif Pinelis to this last post shows that my conjecture was trivially wrong. Anyway, as suggested by Benoît Kloeckner in his insightful comments, this was due to the fact I had assumed random uniform reflection on $z=0$ in my two previous posts. Replacing the random uniform reflection with the Knudsen law should finally give the right answer. 


NOTE2. Let me add here a final remark about the physical meaning of the assumption of Knudsen law in our model. It simply says the every little region of the face $z=0$ emits light with constant radiance in any direction (such an ideal isotropic emitter is often called a Lambertian emitter: see [Lambert's cosine law](https://en.wikipedia.org/wiki/Lambert%27s_cosine_law). Such a behavior is exactly what we need in the thermodynamical situation (blackbody radiation in cavity in thermodynamical equilibrium) considered in the proof of Wien's displacement law (see Planck, The Theory of Heat Radiation, Part II, Section III, especially $\S 78$).



A Really Simple Stochastic Dynamic Billiard
This answer is a modification of my answer at [2D billiard](https://mathoverflow.net/questions/285783/a-really-simple-stochastic-dynamic-billiard). 


To avoid trivialities, assume that the initial position of the particle is in the interior of the box. 
Also, at this point, assume that the first time the particle hits the side $z=0$, it occurs not at a vertex of the box. 
If $x$- and $y$-coordinates of the initial velocity are $0$, then the limit in question is obviously $0$. 


Otherwise, in a nonrandom finite time, depending only on the initial position and velocity, the particle will hit the side $z=0$. So, without loss of generality let us assume the movement begins right after that hitting time, with the initial velocity of magnitude $v$ at a random angle $\theta\_1$ with the $z$-axis. Note that the vertical coordinate of the velocity may change in magnitude only right after hitting the side $z=0$; except for this, it may only change in sign right after hitting the side $z=c$. 
Also, the cosine of the angle between the velocity and the $z$-axis does not vary while the particle is traveling from the side $z=0$ toward the side $z=c$. 


Introduce $Y\_j:=\cos\theta\_j$, $X\_i:=\frac{c}{vY\_i}$ -- the $i$th traveling time from the side $z=0$ to the side $z=c$ (which is equal to the $i$th traveling time from the side $z=c$ back to the side $z=0$), $T\_j:=\sum\_1^j X\_i$, 
$N\_t:=\sum\_{j=1}^\infty I\{T\_j\le t\}$ (with $I\{\cdot\}$ denoting the indicator function). Do not confuse $N\_t$ with $N(t)=\sum\_{j=1}^\infty I\{2T\_{j-1}+X\_j\le t\}$ -- the number of times the particle hits the side $z=c$. 
Let then
\begin{gather\*}
R\_t:=\sum\_{j=1}^{N\_t}Y\_j, \tag{1} \\ 
 S\_t:=\sum\_{j=1}^{N(t)} \cos \theta\_j=\sum\_{j=1}^{N(t)}Y\_j \tag{2} 
 =\sum\_{j=1}^\infty Y\_j\,I\{2T\_{j-1}+X\_j\le t\}, \\
 S\_t^-:=\sum\_{j=1}^\infty Y\_j\,I\{2T\_j\le t\}=R\_{t/2}, \quad
 S\_t^+:=\sum\_{j=1}^\infty Y\_j\,I\{2T\_{j-1}\le t\}. 
\end{gather\*}
Then  

\begin{equation\*}
 S\_t^-\le S\_t\le S\_t^+; \tag{3} 
\end{equation\*}
everywhere here $t$ is a positive real number. 
By the [Renewal Reward Theorem](http://www.randomservices.org/random/renewal/Reward.html), we have 
\begin{equation\*}
 \frac1t\,ER\_t\underset{t\to\infty}\longrightarrow\frac{EY\_1}{EX\_1}=\frac{2/3}{2c/v}=\frac{v}{3c}. 
\end{equation\*}
So, 
\begin{equation\*}
 \frac1t\,ES\_t^-=
 \frac1t\,ER\_{t/2}\underset{t\to\infty}\longrightarrow\frac{v}{6c}. \tag{4} 
\end{equation\*}
Next, by the independence of $Y\_j$ and $T\_{j-1}$ and, again, by the Renewal Reward Theorem (with $1$ in place of $Y\_j$), 
\begin{multline\*}
ES\_t^+=\sum\_{j=1}^\infty EY\_j\,EI\{2T\_{j-1}\le t\}
=EY\_1\,E\sum\_{j=1}^\infty I\{2T\_{j-1}\le t\}
=EY\_1\,(1+EN\_{t/2}) \\ 
= EY\_1\,\Big(1+E\sum\_{j=1}^{N\_{t/2}}1\Big) 
= EY\_1\,\Big(1+\frac t{2+o(1)}\, \frac{E1}{EX\_1}\Big)
\sim EY\_1\,\frac t2\, \frac1{EX\_1}
=\frac{tv}{6c} \tag{5} 
\end{multline\*} 
as $t\to\infty$. 
Thus, by (2)--(5), 
\begin{equation\*}
 \frac1t\,E\sum\_{j=1}^{N(t)} \cos \theta\_j=\frac1t\,ES\_t\underset{t\to\infty}\longrightarrow\frac{v}{6c}. 
\end{equation\*}


