Merging MODIS data using R to produce either single raster layer or multiple raster layers displaying NPP over time

Question

I am trying to collate NPP MODIS data (MOD17A3HGF V061) for an area over Europe (at least 17 tiles) and spanning for the available timeframe (2001-present) in R. I have downloaded the data as .hdf files. The years are in the filenames (see bottom of question for example filename) and when I have read them into ArcPro it is able to recognise that, but it can't mosaic them for some reason. I have tried using various R packages such as MODIStsp and MODISTools but my company's security restricts third party applications downloading data (hence why I have completed the download manually). From what I have been able to achieve (or not as the case is) these tools don't read in the .hdf files, they just work from the online repository (please correct me if wrong). And so I've not been able to use them.

I have thus tried various packages which read .hdf files such as ncdf4, rhdf5, RNetCDF, raster, and terra. When using any of the first four packages, I have read in the list of file names as follows:

library(raster)
library(rhdf5)
library(ncdf4)
library(RNetCDF)

input_directory = "path/to/folder/containing/modis/files/over/several/tiles/and/years"
data_list = list.files(input_directory, pattern = ".hdf", full.names = TRUE, include.dirs = TRUE)

hdfs_rasters <- lapply(data_list, raster) #appeared to work but only a list of file names
# OR
hdfs_nc_open <- lapply(data_list, nc_open) #appeared to work but only a list of file names
# OR
hdfs_hdfile <- lapply(data_list, H5Acreate) # Tis here I learnt the directory is hdf4 files
# OR
hdfs_netcdf <- lapply(data_list, open.nc) #produced error; too many files opened at once

I finally used terra::rast which collected data_list into a list of SpatRasters, and then I created a SpatRasterCollection which I was then able to merge (see code below):

library(terra)

hdfs_rasters <- lapply(data_list, rast)

hdfs_collection <- sprc(hdfs_rasters)

hdfs_merge <- merge(hdfs_collection)

This has merged the tiles spatially (yey!) but I can't extract the temporal element.

I was wondering if there was a way to produce a single raster with the temporal element (NPP/year) as individual bands or if I could produce multiple rasters for the years?

example file name: MOD17A3HGF.A2022001.h19v02.061.2023035051843.hdf

MOD17A3HGF data can be downloaded from here:Earthdata Search on NASA.gov website