It isn't quite false advertising. 

I will attribute it to an omission based on the needs of the many rather than on the needs of the few.

The story is in code...

# ---- start with an array ------

a  # ---- any data, floats for a tiny place over 5 days -----
 
array([[[  0.,   1.,   2.,   3.],
        [  4.,   5.,   6.,   7.],
        [  8.,   9.,  10.,  11.]],

       [[ 12.,  13.,  14.,  15.],
        [ 16.,  17.,  18.,  19.],
        [ 20.,  21.,  22.,  23.]],

       [[ 24.,  25.,  26.,  27.],
        [ 28.,  29.,  30.,  31.],
        [ 32.,  33.,  34.,  35.]],

       [[ 36.,  37.,  38.,  39.],
        [ 40.,  41.,  42.,  43.],
        [ 44.,  45.,  46.,  47.]],

       [[ 48.,  49.,  50.,  51.],
        [ 52.,  53.,  54.,  55.],
        [ 56.,  57.,  58.,  59.]]])

# ---- save to a raster *.tif ----------------

r = arcpy.NumPyArrayToRaster(a)
r.save("c:/temp/np2rast.tif")

# ---- load it back in to check your work ----

b = arcpy.RasterToNumPyArray("c:/temp/np2rast.tif")

b  # ---- seems to only read the first 'day' ---
array([[  0.,   1.,   2.,   3.],
       [  4.,   5.,   6.,   7.],
       [  8.,   9.,  10.,  11.]])

# ---- well that isn't good, hmmmm  ---------------
# The 
rasters = []
for i in range(a.shape[0]):
    ai = a[i]
    rast = arcpy.NumPyArrayToRaster(ai)
    r_name = "in_memory/a{:0>3}.tif".format(i)
    rasters.append(r_name)
    rast.save(r_name)
rasters = ";".join([i for i in rasters])

rasters
'in_memory/a000.tif;in_memory/a001.tif;in_memory/a002.tif;in_memory/a003.tif;in_memory/a004.tif'

fname = "c:/temp/np2rast_comp.tif"
arcpy.management.CompositeBands(rasters, fname)  # ---- let's give Arc* another chance
c = arcpy.RasterToNumPyArray(fname)

c  # ---- everything is back as it should be -----------
 
array([[[  0.,   1.,   2.,   3.],
        [  4.,   5.,   6.,   7.],
        [  8.,   9.,  10.,  11.]],

       [[ 12.,  13.,  14.,  15.],
        [ 16.,  17.,  18.,  19.],
        [ 20.,  21.,  22.,  23.]],

       [[ 24.,  25.,  26.,  27.],
        [ 28.,  29.,  30.,  31.],
        [ 32.,  33.,  34.,  35.]],

       [[ 36.,  37.,  38.,  39.],
        [ 40.,  41.,  42.,  43.],
        [ 44.,  45.,  46.,  47.]],

       [[ 48.,  49.,  50.,  51.],
        [ 52.,  53.,  54.,  55.],
        [ 56.,  57.,  58.,  59.]]])

Now if you followed so far, I have made reference to working with numpy arrays to do a lot of work.  Rarely do I have the need to bring a 'raster' (ie *.tif) back in to ArcGIS PRO or ArcMap.  I did and discovered, that arcpy.RasterToNumPyArray seems to only like bringing back  a limited amount of data, probably with the expectation that user wants to map floating point values (ie 1 band) or integer values ( ie 3 bands for a 'picture'/'map').

This is a work-around, built-in.

You can of course store anything numeric in a tiff.  How about some statistics for our array 'a'.

mins = np.min(a, axis=0)
maxs = np.max(a, axis=0)
avg = np.mean(a, axis=0)
rang = np.ptp(a, axis=0)
medn = np.median(a, axis=0)

stats = [mins, maxs, avg, rang, medn]
rasters = []
for i in range(5):
    ai = stats[i]
    rast = arcpy.NumPyArrayToRaster(ai)
    r_name = "in_memory/a{:0>3}.tif".format(i)
    rasters.append(r_name)
    rast.save(r_name)
rasters = ";".join([i for i in rasters])

arcpy.management.CompositeBands(rasters, "c:/temp/stats.tif")
<Result 'c:\\Temp\\stats.tif'>

s = arcpy.RasterToNumPyArray("c:/temp/stats.tif")

s  # ---- in order, min, max, avg, range, median ----
array([[[  0.,   1.,   2.,   3.],
        [  4.,   5.,   6.,   7.],
        [  8.,   9.,  10.,  11.]],

       [[ 48.,  49.,  50.,  51.],
        [ 52.,  53.,  54.,  55.],
        [ 56.,  57.,  58.,  59.]],

       [[ 24.,  25.,  26.,  27.],
        [ 28.,  29.,  30.,  31.],
        [ 32.,  33.,  34.,  35.]],

       [[ 48.,  48.,  48.,  48.],
        [ 48.,  48.,  48.,  48.],
        [ 48.,  48.,  48.,  48.]],

       [[ 24.,  25.,  26.,  27.],
        [ 28.,  29.,  30.,  31.],
        [ 32.,  33.,  34.,  35.]]])

# ---- It won't make for a nice map, but your data is there

Next up, explorations of the limitations of working with *.tif files for multidimensional raster/array data.  

There are a number of options within the Esri suite of software, some of it relies in part on open source offerings.  Some of these (ie libtiff and tifffile) have even made it into open source software too (ie GDAL, scipy.ndimage and skimage etc etc), kindof like open source-open source.

Enough for now.

This comes up time and again. ....I have 'X' rasters for an area and I need to determine Y over time....

Here is an example.... ....calculating percentiles....

Cell Statistics... is the obvious choice.  But you don't have the Spatial Analyst extension? Doesn't matter.  All the tools that you need are already provided if you are using current versions of ArcMap or ArcGIS PRO.

Do some work

(1)  Reading the raster files

import numpy as np
import arcpy

arrs = []
folder = r"C:\Data\rasters"  # ---- specify a folder containing all the rasters ----
arcpy.env.workspace = folder
rasters = arcpy.ListRasters("*", "TIF")
for raster in rasters:
    arrs.append(arcpy.RasterToNumPyArray(raster))

a = np.array(arrs) # ---- the master array ----

Once the rasters are all converted to arrays, they were appended to a list.  A master array was created in line 9.

As an example, 31 rasters in tif format were read in.  These are small rasters, but they can be as large as you can process.  You are now ready to calculate values for them using numpy.

(2)  Getting the statistics

median = np.median(a, axis=0) # ---- pick a statistic ----
 
array([[ 59.,  54.,  36., ...,  57.,  46.,  46.],
       [ 43.,  45.,  59., ...,  38.,  51.,  51.],
       [ 35.,  50.,  57., ...,  47.,  55.,  65.],
       ..., 
       [ 43.,  52.,  40., ...,  56.,  62.,  60.],
       [ 45.,  57.,  39., ...,  45.,  44.,  48.],
       [ 49.,  57.,  50., ...,  56.,  50.,  58.]]

mins = np.min(a, axis=0)

array([[1, 3, 0, ..., 6, 4, 5],
       [1, 5, 6, ..., 6, 7, 2],
       [3, 2, 4, ..., 2, 5, 4],
       ..., 
       [4, 0, 1, ..., 3, 4, 6],
       [0, 0, 0, ..., 5, 1, 0],
       [4, 1, 1, ..., 0, 3, 7]])

# ---- name a statistic, with or without nodata values, it can be done ----

Now when you are done your calculations you will probably want to make a raster so you can bask in the beauty of your efforts.

(3)  Reading the raster information and saving the result

rast = arcpy.Raster(r01)  # ---- read the first raster we loaded ----

rast = arcpy.Raster(r01)  # ---- simple 

dir(rast)                 # ---- get raster information... some snipping here ----

[..., 'bandCount', 'catalogPath', 'compressionType', 'extent', 'format', 'hasRAT',
 'height', 'isInteger', 'isTemporary', 'maximum', 'mean', 'meanCellHeight',
 'meanCellWidth', 'minimum', 'name', 'noDataValue', 'path', 'pixelType', 'save',
 'spatialReference', 'standardDeviation', 'uncompressedSize', 'width']

# ---- Save the result out to a new raster ------

r01 = rasters[1]                    # --- rasters have the same extent and cell size
rast = arcpy.Raster(r01)

lower_left = rast.extent.lowerLeft  # ---- needed to produce output
cell_size = rast.meanCellHeight     # ---- we will use this for x and y

f = r"c:\temp\result.tif"
r = arcpy.NumPyArrayToRaster(a, lower_left, cell_size, cell_size)
r.save(f)

(4) Saving and reloading the arrays

Now, the nice thing about arrays is that you can save a 3D to a file for reloading later (or any dimension for that matter).

If you have a 3D array like we have, this is kind of like a big 'zip'.  

np.save(r"c:\temp\myarray.npy", a)  # ---- 'a' is the multi dimensional array

a.shape  # ---- 31 days, 100 x 100 cells 
(31, 100, 100)

# ---- reload it 

im_back = np.load(r"c:\temp\myarray.npy")

im_back.shape
(31, 100, 100)

# ---- check to see if it is good

np.all(a == im_back)  # ---- some numpy magic ;)
True

That's about all.  Not covered here is nodata cells (easy, you can assign one) and statistics for nodata (easy, every stat has a nodata equivalent).

Try your skills with the attached *.npy file.