Inelastic Particle Clusters from Cumulative Momenta

Kevin Chien, Aidan Mager, Laurel Safranek, Jackson Zariski

Abstract

We consider a physical system comprising discrete massive particles on the real line whose trajectories interact via perfectly inelastic collision. It turns out that polygons formed in a “cumulative momentum diagram” of the initial conditions allow us to easily predict how many particle clusters form as time $t\to\infty$ . We explore an application of this to a unit mass system with $\pm 1$ velocities.

1 Introduction

Inelastic particles, or sticky particles, have been studied in relation to the pressureless Euler system in $(0,\infty)\times\RR^{d}$ :

\begin{cases}\partial_{t}\rho+div(\rho v)=0\\ \partial_{t}(\rho v)+div(\rho v\otimes v)=0\end{cases}

(1.1)

where $\rho(t,x)$ and $v(t,x)$ are the density and velocity fields of gas particles, respectively [7]. Note that these equations may be derived from the usual Euler equations for ideal compressible fluids by letting the pressure go to zero [2]. This system governs how an ideal pressureless system of gas particles evolves over time. Such a system has been proposed by Zeldovic [8] as a simplified model for early formations of galaxies, where individual dust particles floating in space may collide and form newer, bigger particles. Bianchini [7] proved that this system admits a unique sticky particle solution.

In one dimension, solutions to (1.1) can be obtained in the limit of a discrete system of sticky particles (see [4]) often with an interaction potential. In this paper, we focus on such a discrete system with no potential, where particles move in piecewise linear trajectories. Our main result reframes this system in terms of a cumulative momentum diagram, providing a new geometric picture of the situation that simplifies the determination of how many clusters of particles form. As an application, we use this approach to link the problem of determining how often a system of unit mass, unit speed particles on the line form one cluster to a problem concerning simple symmetric random walks. We also provide the code to implement our cumulative momentum approach.

In Section 2, we provide the necessary definitions and main theorem statement. Section 3 contains the proof of our result. Section 4 uses this approach on our toy problem, and Section 5 gives a Python code sample.

Acknowledgements

We would like to thank Jayadev Athreya for his support and bringing this topic to our attention; Rowan Rowlands for the helpful reference to [3]; as well as the Washington Experimental Math Lab for facilitating this collaboration.

2 Preliminaries

2.1 Inelastic particle systems

A particle $a:=(m,x,v)\in\RR_{>0}\times\RR\times\RR$ is an ordered triple whose elements denote mass, initial position, and initial velocity, respectively. A system of $n$ particles denotes an $n$ -tuple of particles $(a_{i})_{i=1}^{n}$ satisfying $x_{i}\neq x_{j}$ for $i\neq j$ , and ordered such that $x_{i}<x_{j}$ if and only if $i<j$ .

To each particle $a_{i}$ in our system, we may associate a unique trajectory $\gamma_{i}:[0,\infty)\to\RR$ , which is a continuous, piecewise linear map satisfying the following properties (see Proposition 2.1 in [4]):

1.

$\gamma_{i}(0)=x_{i},\,\dot{\gamma}(0+)=v_{i}$ ;
2.

If $\gamma_{i}(s)=\gamma_{j}(s)$ , then $\gamma_{i}(t)=\gamma_{j}(t)$ for all $t\geq s$ (stickiness); and
3.

If $\gamma_{i_{1}}(t)=\cdots=\gamma_{i_{\ell}}(t)\neq\gamma_{i}(t)$ at some $t>0$ , then

$\dot{\gamma}_{i_{j}}(t+)=\frac{\sum_{k=1}^{\ell}m_{i_{k}}\dot{\gamma}_{i_{k}}(t-)}{\sum_{k=1}^{\ell}m_{i_{k}}}$

for $i_{j}\in\{i_{1},\ldots,i_{\ell}\}$ (conservation of momentum).

Figure 1: Clusters for unit mass particles with initial velocities

(+1,-1-1,+1,-1,+1)

In the above, we employ the notation $\dot{\gamma}(t\pm):=\lim_{h\to 0^{\pm}}(\gamma(t+h)-\gamma(t))/h$ . We denote $\dot{\gamma}(t+)$ the velocity of $a_{i}$ at time $t$ : in particular, the velocity exists for all $t\in[0,\infty)$ , but may not be continuous.

Observe that the continuity of trajectories and Property 2 together imply that if particles $a_{i},a_{j}$ collide by time $t$ , then all particles $a_{k}$ for $i\leq k\leq j$ must have collided by time $t$ . Also, by expanding the right hand side of the expression in Property 3, we can recover the velocity of a particle at any time from the initial conditions of all particles involved in the collision. That is, at a given time $t\in[0,\infty)$ , if $J\subset(1,\ldots,n)$ is the subset of indices such that $\gamma_{i}(t)=\gamma_{j}(t)$ for all $j\in J$ , then the velocity of $a_{i}$ is

\displaystyle\dot{\gamma}_{i}(t+)=\frac{\sum_{j\in J}m_{j}v_{j}}{\sum_{j\in J}m_{j}}.

(2.1)

Using these trajectories to define an equivalence relation, we can partition our system into clusters: two particles $a_{i},a_{j}$ are in the same cluster $A_{k}$ if and only if $\lim_{t\to\infty}|\gamma_{i}(t)-\gamma_{j}(t)|=0$ . By our observation above, we may order our clusters $(A_{i})_{i=1}^{\ell}$ such that $i<j$ if and only if the initial position of any particle in $A_{i}$ is less than the initial position of any particle in $A_{j}$ (see Figure 1). The velocity of a cluster is the velocity of any constituent particle, for $t$ large enough such that all particles in the cluster share the same trajectory.

Figure 2: The cumulative momentum diagram for the same system.

2.2 Cumulative momentum diagrams and polygons

Given a system $(a_{i})_{i=1}^{n}$ , we now give an associated construction using initial momenta as follows. Let $P_{k}:=(\sum_{i=1}^{k}m_{i},\sum_{i=1}^{k}m_{i}v_{i})$ , $P_{0}:=(0,0)$ , and let $L_{k}$ be the line segment in $\RR^{2}$ connecting $P_{k-1},P_{k}$ , for $1\leq k\leq n$ . We may view $\Gamma(f):=\bigcup_{k=1}^{n}L_{k}$ as the graph of a continuous, piecewise linear function $f:[0,\sum_{i=1}^{n}m_{i}]\to\RR$ , $f(0)=0$ .

Recall that the convex envelope of $f$ is defined by

conv(f):=\sup\{g\leq f:\,g\,\text{ is convex}\}.

The properties of $f$ imply $r:=conv(f):[0,\sum_{i=1}^{n}m_{i}]\to\RR$ is also a continuous, piecewise linear function; let $\Gamma(r)$ denote its graph and $(R_{i})$ its linear components of constant slope, ordered by increasing $x$ -coordinate of its left endpoint. Observe that $(R_{i})$ has monotonically increasing slope due to the convexity of $r$ . Together, we refer to $\Gamma(f)\cup\Gamma(r)$ as the cumulative momentum diagram of our system.

For each $R_{i}$ , there are points $(P_{i_{k}})$ , with $P_{i_{k}}\in R_{i}$ , ordered by increasing $x$ -coordinate. A polygon $T$ in our diagram refers to either

(a)

The simple polygon bounded by $\Gamma(f)$ and the portion of $R_{i}$ joining $P_{i_{k}},P_{i_{k+1}}$ , for $i_{k+1}\neq i_{k}+1$ ; or
(b)

The line segment $L_{i_{k+1}}$ when $i_{k+1}=i_{k}+1$ (the “degenerate” case).

We can now decompose our diagram into a sequence of polygons $(T_{i})_{i=1}^{\ell}$ , ordered using the $x$ -coordinate of the vertex (or endpoint) with smallest $x$ -coordinate in each polygon. This further induces a partition on our system $(a_{i})_{i=1}^{n}$ , where $a_{i},a_{j}\in T_{k}$ if and only if $L_{i},L_{j}$ belong to an edge of $T_{k}$ . In an abuse of language, the slope of $T$ refers to the slope of the linear component $R$ that bounds $T$ . See Figure 2 for an example diagram, and note $T_{3}$ is a line segment that counts as a polygon by our rules above.

2.3 Statement of theorem

Our main result states that given a system $(a_{i})_{i=1}^{n}$ , there is a bijection between ordered clusters $A_{i}$ and ordered polygons, with $A_{i}=T_{i}$ .

Theorem.

Let $(a_{i})_{i=1}^{n}$ be a system whose cumulative momentum diagram decomposes into polygons $(T_{i})_{i=1}^{\ell}$ . Then the system partitions into clusters $(A_{i})_{i=1}^{\ell}$ with $A_{j}=T_{j}$ , and the velocity of $A_{j}$ is given by the slope of $T_{j}$ .

3 Proof of main result

The key step in showing our main result is the ability to consider polygons in isolation, where particles outside a given polygon never interact with particles within. This is the statement of the following “cutting” lemma, which utilizes the monotonicity of $(R_{i})$ :

Lemma 3.1.

Particles in $T_{i}$ never collide with particles in $T_{\ell}$ for $i\neq\ell$ .

Proof.

By our observation following Property 2, it suffices to show that $\dot{\gamma}_{j}(t+)\leq\dot{\gamma}_{j+1}(t+)$ for all $t\geq 0$ , where $j$ is the index of the particle $a_{j}\in T_{i}$ with largest $x$ -coordinate. Observe that until $a_{j},a_{j+1}$ collide, their respective velocities at any time are given by

	$\displaystyle\dot{\gamma}_{j}(t+)$	$\displaystyle=\frac{\sum_{k=0}^{n_{1}-1}m_{j-k}v_{j-k}}{\sum_{k=0}^{n_{1}-1}m_{j-k}}$		(3.1)
	$\displaystyle\dot{\gamma}_{j+1}(t+)$	$\displaystyle=\frac{\sum_{k=1}^{n_{2}}m_{j+k}v_{j+k}}{\sum_{k=0}^{n_{2}}m_{j+k}}$		(3.2)

for $n_{1},n_{2}\geq 1$ , according to (2.1). The right hand side of (3.1) is the slope of the line passing through $P_{j-n_{1}}$ and $P_{j}$ , which is sharply bounded above by $R$ , the slope of $T_{i}$ ; this follows from $P_{j}\in\Gamma(r)$ and the convexity of $r$ . Similarly, the right hand side of (3.2) is the slope of the line passing through $P_{j}$ and $P_{j+n_{2}}$ , which is sharply bounded below by $\tilde{R}$ , the slope of $T_{i+1}$ (see Figure 3). Due to monotonicity of polygon slopes, we know $R\leq\tilde{R}$ , and therefore $\dot{\gamma}_{j}(t+)\leq\dot{\gamma}_{j+1}(t+)$ as needed. ∎

Figure 3:

Proof of Theorem: We proceed by induction. Suppose our system partitions into one polygon $T$ with slope $R=\frac{\sum_{k=1}^{n}m_{k}v_{k}}{\sum_{k=1}^{n}m_{k}}$ . Let $(A_{i})_{i=1}^{N}$ be the clusters of our system, and denote the velocity of $A_{i}$ by $W_{i}$ . If $J$ is the maximal collection of indices of particles belonging to $A_{i}$ , let $M_{i}:=\sum_{k\in J}m_{k}$ .

Because these clusters do not collide, we have $W_{i}\leq W_{i+1}$ . The velocity $W_{1}$ is given by (3.2) with $j=0$ , which by the reasoning in the proof of Lemma 3.1 gives $W_{1}\geq R$ . Combined with conservation of momentum, this yields

\sum_{k=1}^{n}m_{k}v_{k}=\sum_{j=1}^{N}M_{j}W_{j}\geq W_{1}\sum_{j=1}^{N}M_{j}\geq R\sum_{k=1}^{n}m_{k}=\sum_{k=1}^{n}m_{k}v_{k}.

It follows immediately that $W_{1}=R$ and $N=1$ .

Assume our theorem holds for $\ell-1$ polygons, and consider a system partitioned into $\ell$ polygons. Let $a_{j}$ be the particle with largest $x$ -coordinate belonging to the $T_{\ell-1}$ . By Lemma 3.1, particles in $\bigcup_{i=1}^{\ell-1}T_{i}$ do not collide with particles in $T_{\ell}$ . Thus, the first $\ell-1$ polygons may be considered in isolation from the rest of the system, so by our inductive hypothesis, $(a_{i})_{i=1}^{j}$ partitions into $\ell-1$ clusters with $A_{j}=T_{j}$ . We are then left with the particles in $T_{\ell}$ , which upon applying our base case reasoning corresponds to the single cluster $A_{\ell}$ . ∎

4 Application to a toy problem

Consider the following problem: Given a system of $n$ particles of unit mass with velocities taking values in $\{\pm 1\}$ , we want to determine the probability of one cluster forming. By looking at the associated cumulative momentum diagrams, this is equivalent to summing over the probabilities that a simple symmetric random walk (SSRW) on $\mathbb{Z}$ starting at $0$ and ending at $c\in\mathcal{V}$ stays strictly above $\frac{c}{n}j$ for each step $0<j<n$ , where $\mathcal{V}:=\{-n\leq c\leq n:c=n-2k,k\in\ZZ\}$ is the range of possible total velocities of the $n$ particles.

This type of problem appears in combinatorics in relation to the classic ballot problem and counting North-East lattice paths – or staircase walks – on $\ZZ^{2}$ (see Chapters 1, 2 in [6]). To see this, observe that each SSRW satisfying our needs above corresponds to a lattice path on $\ZZ^{2}$ starting at the origin and ending at $(x,y)=((n-c)/2,(n+c)/2)$ that stays strictly above the line of slope $y=x((n+c)/(n-c))$ ; here, each unit step in the $+x$ (resp. $+y$ )-direction corresponds to a $-1$ (resp. $+1$ ) step in our SSRW.

4.1 Young diagrams and Narayana’s path counting

In order to approach the situation above, we briefly recall Young diagrams and a useful path counting formula. Given $\lambda=(\lambda_{1},\ldots,\lambda_{k})\in\ZZ_{+}^{k}$ with $\lambda_{i}\geq\lambda_{i+1}$ the Young diagram of shape $\lambda$ comprises left justified rows with row $i$ consisting of $\lambda_{i}$ unit boxes.

Let us introduce the notation

{\binom{n}{k}}_{+}:=\begin{cases}{\binom{n}{k}}&\text{\, if\, }k\geq 0;\\ 0&\text{\, if\, }k<0.\\ \end{cases}

Figure 4: Example of a North-East lattice path on the Young diagram of shape

\lambda=(4,2,1,0,0)

We know the following formula:

Theorem (Path Counting, [3] [5]).

The number $N(\lambda)$ of North-East lattice paths joining the southwest corner to the northeast corner of a Young diagram of shape $\lambda$ is given by

\displaystyle N(\lambda)=\det\left({\binom{\lambda_{j}+1}{j-i+1}}_{+}\right)_{1\leq i,j\leq k}

(4.1)

Observe that the matrix that we are taking the determinant of in (4.1) has all ones on the main subdiagonal and zeros all below that. Also, if we extend the notion of Young diagram to tuples $\lambda\in\ZZ_{\geq 0}^{k}$ with the zero entries corresponding to rows with left-aligned single points instead of boxes, equation (4.1) still holds due to determinant properties of block lower-triangular matrices; in this case, $N(\lambda)$ counts the paths on the diagram given by the corresponding Young diagram for nonzero entries of $\lambda$ with the additional vertical line segment joining the single points (Figure 4).

4.2 Solution to the problem

Observe that for $c=\pm n,\pm(n-2)$ , the probabilities of one cluster forming are $0$ and $(1/2)^{n}$ respectively. For $c\in\mathcal{V}$ that are not those values, the number of $\ZZ^{2}$ lattice paths from the origin to $((n-c)/2,(n+c)/2)$ staying above $y=x((n+c)/(n-c))$ is given by $N(\lambda)$ for

\lambda=\left(\left\lceil\left(\frac{n-c}{n+c}\right)\left(\frac{n+c}{2}-j\right)-1\right\rceil\right)_{j=1}^{\frac{n+c}{2}-1}.

This $\lambda$ is the shape of the maximal Young diagram that fits inside $[0,(n-c)/2]\times[0,(n+c)/2]$ and stays strictly above $y=x((n+c)/(n-c))$ . Plugging in $c=-n+2k$ for $2\leq k\leq n-2$ , we find the probability of one cluster forming to be

(1/2)^{n-1}+\frac{1}{2^{n}}\sum_{k=1}^{n-1}\det\left({\binom{\left\lceil\left(\frac{n-k}{k}\right)\left(k-j\right)-1\right\rceil+1}{j-i+1}}_{+}\right)_{1\leq i,j\leq k-1}.

Noting in this case that the number of permissible lattice paths from the origin to $(p,q)$ is the same as that for $(q,p)$ (see Figure 5), we may rewrite the above slightly and conclude that the probability of one cluster forming for $n>2$ even (resp. odd) is:

	$\displaystyle\frac{1}{2^{n-1}}+\frac{1}{2^{n}}\det\left({\binom{\left\lceil(n/2)-j-1\right\rceil+1}{j-i+1}}_{+}\right)_{1\leq i,j\leq(n/2)-1}$
	$\displaystyle+\frac{1}{2^{n-1}}\sum_{k=2}^{n/2-1}\det\left({\binom{\left\lceil\left(\frac{n-k}{k}\right)\left(k-j\right)-1\right\rceil+1}{j-i+1}}_{+}\right)_{1\leq i,j\leq k-1}$

\left[resp.\,\hskip 5.0pt\frac{1}{2^{n-1}}+\frac{1}{2^{n-1}}\sum_{k=2}^{\frac{n-1}{2}-1}\det\left({\binom{\left\lceil\left(\frac{n-k}{k}\right)\left(k-j\right)-1\right\rceil+1}{j-i+1}}_{+}\right)_{1\leq i,j\leq k-1}\right]

Figure 5:

4.3 Some remarks on asymptotics

We conclude this section with a brief discussion on the asymptotics of the probability that one cluster forms from our unit mass and speed system as the number of particles $n\to\infty$ . Our approach to find the above formula does not appear to be conducive in extracting asymptotics due to the complicated determinant terms arising from path-counting. Some relevant asympototics for path-counting below a certain class of lines of rational slope can be found in [1], though the techniques there do not seem to easily generalize.

Note that for a SSRW $(S_{k})_{k=1}^{2n}$ starting at $0$ , the probability $S_{k}$ remains strictly positive at each step given the walk ends at $0$ is

\PP\left(\bigcap_{1<k<2n}(S_{k}>0)|S_{2n}=0\right)=\left(\frac{1}{2}\right)^{2n}C_{n-1}=\left(\frac{1}{2}\right)^{n}\frac{1}{n}\binom{2n-2}{n-1}\sim\frac{e}{2\sqrt{\pi}}\left(\frac{1}{\sqrt{n}}\right),

where $C_{k}$ denotes the $k$ -th Catalan number. Applying Stirling’s formula once again to find

\PP(S_{2n}=0)\sim\frac{1}{\sqrt{\pi}}\left(\frac{1}{\sqrt{n}}\right),

we get

\PP\left(\bigcap_{1<k<2n}(S_{k}>0),S_{2n}=0\right)\sim\frac{e}{2\pi}\left(\frac{1}{n}\right).

Heuristically, because we expect most of our walks to satisfy $S_{2n}=c$ for small values of $c$ with high probability, we conjecture that summing over these possibilities leads to a similar asymptotic of $\sim C/n$ .

5 Computational applications

We seek an algorithm that takes in a sequence of original particles as a list of $(M_{i},V_{i})$ pairs, and outputs a similar list of clusters after all collisions occur. An outline for one implementation of this code utilizing recursion is as follows:

1.

Begin with an initial cumulative momentum diagram for the entire sequence. For each point on the diagram in which the slope between it and its previous point is less than or equal to the slope between it and the following point, add said point to the list. End the list with the final point of the diagram.
2.

Recursively repeat the previous step with just these points acting as a revised, shortened momentum diagram. The terminating condition for said recursion occurs when this process no longer produces a shortened diagram. Take this constant result as the final diagram.
3.

The line segments between consecutive points in this final diagram represent the final polygons in the diagram–the slopes of the segments are the velocities of the final particles while the horizontal magnitude of the line is the mass in the unit case.

This algorithm results in a method of determining the final state of coalescence of a sequence in a much more streamlined manner than the brute-force method involving purely the physical calculations. This proves to be especially helpful in the construction of probability density functions for the state of an $n$ -length system as time goes to infinity. In this type of problem, we seek to find the probabilities of there being $k\leq n$ final clusters in an $n$ -length unit-mass system when assigned random velocities of $\pm 1$ . For low values of $n$ one can fairly easily compute these functions by hand, as each $n$ -length sequence encompasses a total of $2^{n}$ possible configurations. However, to analyze enough of the lower $n$ -length sequences, this exponential growth in the size of the input requires an algorithm that can save as much time as possible when compared to that using purely physical calculations. Following trials on randomly generated sets of sequences, algorithms based on analyzing momentum diagrams with the polygon theorem as as a basis achieves this goal; and, in addition, run with more economical memory usage as well. An example of Python code utilizing this algorithm to find the final state of an arbitrary sequence of particles is below.


# Takes two [x,y] points as input
# Returns slope between both points
def slope(first, second):
    rise = second[1] - first[1]
    run = second[0] - first[0]
    return rise/run

# Takes a list of all [x,y] points of a momentum diagram as input
# Returns only the points used as nodes in said momentum diagram
def polygonAnalysis(seq):
    if len(seq) < 2:
        return seq

    else:
        mins = [seq[0]]
        curr= seq[1]
        currIndex = 1
        while currIndex < len(seq) - 1:
            prev = seq[currIndex - 1]
            next = seq[currIndex + 1]
            local_Min = False
            if slope(prev, curr) <= slope(curr, next):
                local_Min = True
            if local_Min:
                mins.append(seq[currIndex])
            currIndex = currIndex + 1
            curr = seq[currIndex]


        mins.append(seq[len(seq) - 1])
        # Algorithm terminates when no further reductions occur
        if mins == seq:
            return mins
        else:
            return polygonAnalysis(mins)

# Takes an original sequence of particles in the form [mass, velocity] as input
# Returns the final list of particles after all collisions ordered by position.
def findFinalTree(origSequence):
    origParticles = []
    index = 0
    for item in origSequence:
        origParticles.append(item)
        index = index + 1

    numStart = 0
    numDown = 0

    curr_x = numStart
    curr_y = numDown

    rldList = [[curr_x, curr_y]]
    for item in origParticles:
        curr_x = curr_x + item[0]
        curr_y = curr_y + (item[1] * item[0])
        rldList.append([curr_x, curr_y])

    middleParticles = polygonAnalysis(rldList) # Recursive step
    currMIndex = 0
    finalParticles = []
    while currMIndex < len(middleParticles) - 1:
        firstParticle = middleParticles[currMIndex]
        secondParticle = middleParticles[currMIndex + 1]
        mass = secondParticle[0] - firstParticle[0]
        finalParticles.append([mass, slope(firstParticle, secondParticle)])
        currMIndex = currMIndex + 1

    return [finalParticles]


# Runs program, prints to terminal list of final particles after all collisions
origSequence = ANY_ARBITRARY_PARTICLE_LIST
print(findFinalTree(origSequence))

References

[1] Cyril Banderier and Michael Wallner “The kernel method for lattice paths below a line of rational slope” In Lattice path combinatorics and applications Springer, 2019, pp. 119–154
[2] Yann Brenier and Emmanuel Grenier “Sticky Particles and Scalar Conservation Laws” In SIAM J. Numer. Anal. 35.6 USA: Society for IndustrialApplied Mathematics, 1998, pp. 2317–2328 DOI: 10.1137/S0036142997317353
[3] Mihai Ciucu “A short conceptual proof of Narayana’s path-counting formula”, 2016 arXiv:1602.02085 [math.CO]
[4] Ryan Hynd “Sticky particles and the pressureless Euler equations in one spatial dimension”, 2020 arXiv:1804.09340 [math.AP]
[5] Germain Kreweras “Sur une classe de problèmes de dénombrement liés au treillis des partitions des entiers” In Cahiers du Bureau universitaire de recherche opérationnelle Série Recherche 6 Institut Henri Poincaré - Institut de Statistique de l’Université de Paris, 1965, pp. 9–107 URL: http://www.numdam.org/item/BURO_1965__6__9_0/
[6] T.. Narayana “Lattice Path Combinatorics with Statistical Applications; Mathematical Expositions 23” University of Toronto Press, 1979 URL: http://www.jstor.org/stable/10.3138/j.ctvfrxdw5
[7] Bianchini Stefano and Daneri Sara “On the sticky particle solutions to the multi-dimensional pressureless Euler equations”, 2020 arXiv:2004.06557 [math.AP]
[8] Ya B Zel’Dovich “Gravitational instability: An approximate theory for large density perturbations.” In Astronomy and astrophysics 5, 1970, pp. 84–89