You probably used the stock
arima function to model ARMA process. The stock function uses MLE with normal distribution. If you think your errors are bimodal and the process is still ARMA, then the only solution is to get into nonnormal modeling. You may need to implement your own MLE with your own error distribution.
However, looking at your Fig.2 I'd be wary of those three spikes in 2006-2008. They're the ones making your distribution look bimodal. Is it a real phenomenon? Could these by outliers? What's worrying me is that they are accompanied by huge drops, a sort of a zig-zag move. This often indicates an error or irregularity in data. For instance, due to some rule or event you missed some observations, then picked them in next period. In this case you'll see a drop, then spike.
Bimodal distributions are possible, but settling on a bimodal distribution simply based on a histogram in levels is a stretch. I'd first study those zig-zags, at the very least, drop those observations and redraw the histogram. Chances are you'll find something unusual about these observations.