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
```
[![random pure density matrix][1]][1]

```
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][miszcazk] (Miszczak 2011).

 [1]: https://i.sstatic.net/fMuRI.png

[miszcazk]: https://arxiv.org/pdf/1102.4598.pdf