As a "toy" application, we consider a planet with an optically thick, Rayleigh-scattering atmosphere, illuminated by a star (assumed far enough from the planet that we may consider all the incident rays to be paralel).

The radiation incident upon the atmosphere will span the full range of mu_i = cos(theta_i) values, where theta_i is the angle between the incoming radiation and the normal to the surface of the atmosphere. We have thus computed an array of solutions for 22 values of mu_i : mu_i = 0.0174 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.965 0.9986 1. This allows interpolation to any desired mu_i. The total optical depth of the slab was set to 2T=9 (mid-plane depth T=4.5) and the albedo set to unity. At each scattering, the "photon" is split, with the fraction that escapes recorded, while the fraction that does not reach the boundary is just rescattered. Since we keep the emerging Stokes parameters for each of the scatterings, we can obtain the solution for any single-scattering albedo by simply multiplying the nth set of Stokes parameters by albedo^n and summing the results.
The solutions were computed with the code Pol_beam.ijs. We ran groups of 600,000 photons, each scattered 60 times; this was repeated 800 times (i.e., a total of 0.48 billion photons for each of the 22 mu_i cases.) The outgoing photons were binned into 30 values of mu and 36 values of phi (over the range of 0-180 deg only due to symmetry).
Here are results PB_tau_4.5 in a form easily read by a J routine: read_beams45.ijs.

Here is the J code to evaluate the scattered radiation at about 1500 points on the illuminated disk: PLANET.ijs.

It should be pointed out that while the radiation which emerges after one scattering must be polarized perpendicular to the plane defined by the star and the observer, the multiply scattered radiation will have its plane of polarization influenced by the local tilt of the surface. Since the overall tilt of the polarization angle is generally small, the angle has been multiplied by a scaling factor for these illustrations. The length of the blue segment indicates the magnitude of the polarization (scale at upper right). The orientation of the segment indicates the local plane of polarization (but the deviation from 0 deg - vertical - has been multiplied by the "Pol. tilt scale" factor).
Polarization at 30 deg phase angle.
Polarization at 45 deg phase angle.
Polarization at 60 deg phase angle.
Polarization at 75 deg phase angle.
Polarization at 90 deg phase angle.
Polarization at 110 deg phase angle.
Polarization at 125 deg phase angle.
Polarization at 135 deg phase angle.
(The J code to plot these figures is pol_pict.ijs.)

Since the polarization tilts of the upper and lower hemispheres are opposite, and thus cancel, the integrated polarization will be exactly perpendicular to the direction of the illumination. Here are the resultant values of intensity and polarization integrated over the visible disk: Flux and Polarization vs. Phase.

We can average the Stokes parameters over just the upper hemisphere, to see how the tilt varies with phase: Tilt vs. Phase.
The tilt of the lower hemisphere is just opposite, so the net polarization will be exactly vertical. (But when the phase exceeds ~165 deg, the tilt passes 45 deg and the net polarization from both hemispheres switches from vertical to horizontal.)

Finally, here is a comparison of the brightness of the planet as a function of phase compared to the "Lambert sphere": Flux vs Lambert Sphere. (A Lambert sphere is a sphere with a surface such that any incoming beam of light is reflected with equal intensity in all directions.)