One approach is to start with a random pure state and then form a density matrix in that state:

randomPureState[n_] := RandomComplex[{-1-I, 1+I}, n] // Normalize singleStateDensityMatrix[state_] := Outer[Times, state, Conjugate[state]] randomPureDensityMatrix[n_] := singleStateDensityMatrix @ randomPureState @ n 

We can also define a test for validity:

test[m_] := <| "Hermitian" -> HermitianMatrixQ[m] , "PosSemiDef" -> PositiveSemidefiniteMatrixQ[m] , "Trace1" -> Tr[m]==1 , "Idempotent" -> AllTrue[Chop[m.m-m], # == 0 &, 2] |> // <| #, "Valid" -> And@@# |> & 

So then:

SeedRandom[1] $m = randomPureDensityMatrix[4] $m // MatrixForm 

test[$m] (* <|"Hermitian"->True,"PosSemiDef"->True,"Trace1"->True,"Idempotent"->True,"Valid"->True|> *) 

AllTrue[Table[randomPureDensityMatrix[4], 100000], test[#]["Valid"] &] (* True *) 

The present definition of randomPureState neglects the extremely tiny possibility of generating an invalid null state. It is left to the reader to add that check if desired.

One approach is to start with a random pure state and then form a density matrix in that state:

randomPureState[n_] := RandomComplex[{-1-I, 1+I}, n] // Normalize singleStateDensityMatrix[state_] := Outer[Times, state, Conjugate[state]] randomPureDensityMatrix[n_] := singleStateDensityMatrix @ randomPureState @ n 

We can also define a test for validity:

test[m_] := <| "Hermitian" -> HermitianMatrixQ[m] , "PosSemiDef" -> PositiveSemidefiniteMatrixQ[m] , "Trace1" -> Tr[m]==1 , "Idempotent" -> AllTrue[Chop[m.m-m], # == 0 &, 2] |> // <| #, "Valid" -> And@@# |> & 

So then:

SeedRandom[1] $m = randomPureDensityMatrix[4] $m // MatrixForm 

test[$m] (* <|"Hermitian"->True,"PosSemiDef"->True,"Trace1"->True,"Idempotent"->True,"Valid"->True|> *) 

AllTrue[Table[randomPureDensityMatrix[4], 100000], test[#]["Valid"] &] (* True *) 

Caveats

The present definition of randomPureState neglects the extremely tiny possibility of generating an invalid null state. It is left to the reader to add that check if desired.

Also, randomPureState is simplistic and does not generate states uniformly across the state space. If uniformity is desired, then one should use more elaborate methods. See, for example, Generating and using truly random quantum states in Mathematica (Miszczak 2011).