Drawing 3D eigenvectors in R

Question

I'm currently attempting to construct a 3D confidence region based on the following quadratic form;

confidence region for mean mu

while I have the ellipsoid part figured out; I can't seem to figure out a way to include the major and minor axis. As per the textbook I'm using, the major axis is length of major axis units in the direction of $e_i$--the ith eigenvector of the estimate of the covariance matrix S. For plotting the ellipsoid, I used the following code that required the rgl package:
plot3d(ellipsoid(mu_1,S_1,(((n-1)*p)/(n-p))*(qf(0.05,df1=p,df2=n-p,lower.tail = FALSE)),segments=100),alpha=0.5,color='blue')

But I believe the function creates its own 3D mesh object within which it produces the ellipse, so I'm having a lot of trouble constructing the 3D eigenvectors on the same 3D plane. In fact, even if I completely abandon the ellipsoid part, I can't seem to find a way to even just draw out the 3D eigenvectors centered at $\bar{x}$. I've loaded packages like plotly and matlib and even tried something as dumb as defining a 3D parametric function of a line that went as follows:

lines <- function(x0,a){ t <- seq(0,1,length=10000) return(matrix(x0+c(a[1]*t,a[2]*t,a[3]*t),nrow=10000,ncol=3,byrow=TRUE)) }

I then went ahead with the plot3d function, to no avail. Any help would be appreciated!