Can Squeezed Vacuum Produce non-zero $langle E rangle$?

I believe that squeed vacuum can be represented in the Fock state basis as:

$|mathrm{SMSV}rangle=frac{1}{sqrt{cosh r}} sum_{n=0}^{infty}left(-e^{i phi} tanh rright)^{n} frac{sqrt{(2 n) !}}{2^{n} n !}|2 nrangle$

Looking at pictures, it seems as though the electric-field representation looks like this:

Which makes intuitive sense, as you have vacuum that becomes squeezed or unsqueezed based on the phase.

A rotated plot of $E(phi)$ this should look like this:

Where the y-axis is the phase, and the x is the observed E-field amplitude.

When I try to plot the wavefunction for these squeezed states, I get a different plot (I’m plotting a contour plot):

The picture to the left is an example contour plot for a coherent state, while the right plot is a plot of squeezed vacuum.
In this case, you can see that the mean photon number appears to change a bit with phase, and there is also an asymmetry in the noise.

Am I correct to assume that squeezed vacuum should not behave like this (that the expected value of E vs phase should also be 0)? Judging by the contour plot alone, I would’ve guessed that I am plotting a squeezed coherent state.

There is a point in the contour plot that exhibits squeezing at zero amplitude, so maybe that is why it’s called squeezed vacuum?

Here’s my Mathematica code if anyone is interested. (I also looked at a coherent state to confirm that it’s not just a mistake in my code.)

SetOptions[Plot, Frame -> True, Axes -> True, 
  LabelStyle -> {FontFamily -> "Arial", FontSize -> 30}, 
  ImageSize -> {200, 200}, Frame -> True, 
  FrameTicks -> {{None, None}, {{0}, None}}, 
  FrameLabel -> {{None, None}, {None, None}}, 
  GridLinesStyle -> LightGray, BaseStyle -> 12];

Energy[n_] := (2 n + 1) [HBar]/2 [Omega];

[Psi][z_, n_] := 
  1/2 1/Sqrt[
     2^n n!] ((m [Omega])/([Pi] [HBar]))^(1/
      4) Exp[-((m [Omega] z^2)/(2 [HBar]))] HermiteH[n, 
    Sqrt[(m [Omega])/[HBar]] z];

m = 1;
[Omega] = 1;
[HBar] = UnitConvert[Quantity[1, "PlanckConstant"], "SIBase"];
[HBar] = QuantityMagnitude[[HBar]];
[HBar] = 1;

squeezedstate[r_, [Phi]_] := 
  1/Cosh[r] Sum[[Sqrt]Factorial[(2 n)]/(
     2^n n!) (-E^(I [Phi]) Tanh[r])^n [Psi][z, n], {n, 0, 30}];

alphastate[[Alpha]_, [Phi]_] := 
 Sum[([Alpha] E^(I [Phi]))^n/[Sqrt](n!) [Psi][z, n], {n, 0, 15}]

f3 = ContourPlot[{Abs[alphastate[1, [Phi]]]^2}, {z, -4, 
    4}, {[Phi], -[Pi], [Pi]}];

f4 = ContourPlot[{Abs[squeezedstate[1.5, [Phi]]]^2}, {z, -4, 
    4}, {[Phi], -[Pi], [Pi]}];


comboGrid= 
 Grid[{{ Labeled[f3,  "[Alpha](E)", Top, LabelStyle -> Large],  
    Labeled[f4,  "[Zeta](E)", Top, LabelStyle -> Large]}}]

SetOptions[Plot, Frame -> True, Axes -> True, 
  LabelStyle -> {FontFamily -> "Arial", FontSize -> 30}, 
  ImageSize -> {200, 200}, Frame -> True, 
  FrameTicks -> {{None, None}, {{0}, None}}, 
  FrameLabel -> {{None, None}, {None, None}}, 
  GridLinesStyle -> LightGray, BaseStyle -> 12];

Energy[n_] := (2 n + 1) [HBar]/2 [Omega];

[Psi][z_, n_] := 
  1/Sqrt[2^n n!] ((m [Omega])/([Pi] [HBar]))^(1/
      4) Exp[-((m [Omega] z^2)/(2 [HBar]))] HermiteH[n, 
    Sqrt[(m [Omega])/[HBar]] z];

m = 1;
[Omega] = 1;
[HBar] = UnitConvert[Quantity[1, "PlanckConstant"], "SIBase"];
[HBar] = QuantityMagnitude[[HBar]];
[HBar] = 1;

squeezedstate[r_, [Phi]_] := 
  1/Cosh[r] Sum[[Sqrt]Factorial[(2 n)]/(
     2^n n!) (-E^(I [Phi]) Tanh[r])^n [Psi][z, 2 n], {n, 0, 30}];

alphastate[[Alpha]_, [Phi]_] := 
 Sum[([Alpha] E^(I [Phi]))^n/[Sqrt](n!) [Psi][z, n], {n, 0, 15}]

f3 = ContourPlot[{Abs[alphastate[1, [Phi]]]^2}, {z, -4, 
    4}, {[Phi], -[Pi], [Pi]}];

f4 = ContourPlot[{Abs[squeezedstate[1, [Phi] + [Pi]]]^2}, {z, -4, 
    4}, {[Phi], -[Pi], [Pi]}];


FancyGrid = 
 Grid[{{ Labeled[f3,  "[Alpha](E)", Top, LabelStyle -> Large],  
    Labeled[f4,  "[Zeta](E)", Top, LabelStyle -> Large]}}]