In RealityKit, I want to model a sky with stars. I created a sky sphere, and attached little spheres as stars to it.
This works, but for a higher number of stars it is very slow, e.g. for 6000 stars, it takes 20 sec to create the sky sphere (on a M1 MacBook).

I SceneKit, it is possible to render vertices as points of a certain radius. This is very fast. Unfortunately, this is not possible in RealityKit.

Below is my slow code. Of course it makes no sense to render each point, which is only visible as a little white disk, with a sphere mesh. What is a reasonable approach to speed up the code?

Maybe it is possible (but I have only very little experience) to render the spheres with Metal.
I would pass the point coordinates and the point sizes to the GPU that computes a texture that could be mapped onto the interior of the sky sphere.

Slow code:

    /// Creates a sphere entity anchored at the origin, and sets up the stars that should be attached to it.
    /// - Parameter radius: The radius of the sphere
    /// - Returns: The sphere anchor entity
    func createSkySphere(radius: Length) -> AnchorEntity {
        // Extract stars from the catalog, or use test stars
        Sky.setupUsedStars(maxNrStars: maxNrStars, testStars: testStars, coordinate: testStars == .nadir ? deviceCoordinate: nil)
        
        // Create stars and add them to the sky sphere
        var material = UnlitMaterial(applyPostProcessToneMap: false) // Create a material for the sky sphere
        material.color = .init(tint: .black.withAlphaComponent(0.9))
        let mesh = MeshResource.generateSphere(radius: Float(radius.m)) // Create the sky sphere mesh
        let skySphereEntity = ModelEntity(mesh: mesh, materials: [material]) // Create the model entity
        createStars(sphereEntity: skySphereEntity, stars: Sky.usedStars, skyRadius: radius)
        
        let anchorEntity = AnchorEntity(world: .zero) // Create an anchor entity
        anchorEntity.addChild(skySphereEntity) // Add the sphere to it
        
        return anchorEntity
    }
    
    /// Creates star entities using a sphere entity and predefined stars.
    /// The stars are in an array of star categoties, where each category contains an array of stars of comparable brightness..
    /// - Parameters:
    ///   - sphereEntity: The sphere entity to that the stars are attached
    ///   - stars: An array of star categoties with an array of stars
    ///   - skyRadius: The radius of the sky sphere
    func createStars(sphereEntity: ModelEntity, stars: [[Star]], skyRadius: Length) {
        // Find min and max magnitude of all stars
        let allStars = Array(stars.joined())
        let magnitudes = allStars.map { $0.mag }
        let minMagnitude = magnitudes.min()! // Max brightness
        let maxMagnitude = magnitudes.max()! // Min brightness
        
        let r = Float(skyRadius.m) // Radius of the sky sphere
        for index in 0 ..< stars.count {
            let category = stars[index]
            for star in category {
                let magnitude = star.mag
                let normalizedMagnitude: Double
                if maxMagnitude != minMagnitude {
                    normalizedMagnitude = (magnitude - minMagnitude) / (maxMagnitude - minMagnitude) // 0 ... 1
                } else {
                    normalizedMagnitude = 0
                }
                let color = UIColor(white: Sky.brightnessFromMagnitude(normalizedMagnitude), alpha: 1.0)

                // Convert (ra, dec) to spherical coordinates, in the sky sphere local coordinate system.
                /// ???? is the angle in the x,y-plane from the positive x axis towards the positive y axis, i.e. counter clockwise. It is the same as ra.
                /// ???? is the angle from the positive z axis towards the x,y-plane. It is 90° - dec.
                let (????, ????) = Sky.spherical????????(ra: star.ra, dec: star.dec)
                
                // Convert spherical coordinates to Cartesian coordinates, in the sky sphere local coordinate system.
                // A point with ???? = 90°, ???? = 0° is on the positive x axis at x = r.
                // Here, A point with dec = 0, ra = 0° is on the positive x axis at x = r.
                
                // Convert spherical coordinates to Cartesian coordinates
                let x = r * Float(sin(????.radians) * cos(????.radians))
                let y = r * Float(cos(????.radians))
                let z = r * Float(sin(????.radians) * sin(????.radians))
                
                // Create a small sphere for the star
                let starSize = makeStarSize(category: index, nrCategories: stars.count, maxStarSize: 0.05, minStarSize: 0.01)
                let starMesh = MeshResource.generateSphere(radius: starSize)
                var starMaterial = UnlitMaterial(applyPostProcessToneMap: false)
                starMaterial.color = .init(tint: color)
                let starEntity = ModelEntity(mesh: starMesh, materials: [starMaterial])
                
                // Position the star on the sphere surface
                starEntity.position = [x, y, z]
                
                // Add the star to the sky sphere entity
                sphereEntity.addChild(starEntity)
            }
        }
    }
    
    /// Computes a star size, i.e. the radius of a sphere that represents a star on the sky, depending on the star category.
    /// Star category 0 corresponds to max brightness and thus to max star size.
    /// - Parameters:
    ///   - category: The star category
    ///   - nrCategories: The number of star categories
    ///   - maxStarSize: The max star size
    ///   - minStarSize: The min star size
    /// - Returns: The star size
    func makeStarSize(category: Int, 
                      nrCategories: Int, 
                      maxStarSize: Float, 
                      minStarSize: Float) -> Float {
        assert(nrCategories >= 0)
        assert(category >= 0 && category < nrCategories)
        assert(maxStarSize >= minStarSize)
        guard nrCategories > 1 else {
            return maxStarSize
        }
        let starSize = maxStarSize - (maxStarSize - minStarSize) * Float(category) / Float(nrCategories - 1)
        return starSize
    }
}